###############################################
#FIGURES D1 - D6
###############################################

library(estimatr)
library(dplyr)
library(ggplot2)
library(patchwork)

####load cleaned data####
data1 <- read.csv("./data/data1.csv")
data2 <- read.csv("./data/data2.csv")
data3 <- read.csv("./data/data3.csv")

#set values for plotting
z_95 <- 1.96
z_90 <- 1.645
factor_levels <- c("Identity Threat\nvs. Control", "Whataboutism\nvs. Identity Threat", "Whataboutism\nvs. Control")
options(ggplot2.discrete.colour = function() scale_color_viridis_d())

####Methodological note####

# In several subgroup analyses the sample size becomes smaller and some observations have unusually strong influence
# on the fitted regression line. In these cases, the default robust variance estimator can become unstable and
# may occasionally return missing (NA) standard errors due to numerical leverage problems. To ensure valid and stable statistical inference, 
# I therefore use leverage-adjusted heteroskedasticity-consistent estimators (HC3) for all heterogeneous effects analyses, which apply an additional correction
# for influential observations. This is solely in order to have valid confidence intervals for all reported results.

####TABLE D1 - HTEs for Explicit Prejudice####

####TABLE D1(A) - Explicit Prejudice by Left/Right Self-Placement####

#Left-wing
data1_left <- subset(data1, data1$lr_1 < 5)
data2_left <- subset(data2, data2$lr_1 < 5)
data3_left <- subset(data3, data3$lr_1 < 5)

#models by treatment type
lm1a <- lm_robust(prej_direct_1 ~ as.factor(trt) + education + income + newjob, data=data1_left, se_type="HC3")
lm2a <- lm_robust(prej_direct_1 ~ as.factor(trt) + lr_1 + natid_covar1_1 + income + newjob, data = data2_left, se_type="HC3")
lm3a <- lm_robust(prej_direct_1 ~ as.factor(trt) + lr_1 + unemployed, data = data3_left, se_type="HC3")

idt <- c("Identity Threat\nvs. Control", lm1a$coefficients[2])
idt[2] <- lapply(idt[2], as.numeric); idt <- as.data.frame(idt)
colnames(idt) <- c("Treatment_Type", "Effect")
idt$se <- lm1a$std.error[2]
idt <- idt %>% 
  mutate(
    cih_95 = Effect + z_95 * se,
    cil_95 = Effect - z_95 * se,
    cih_90 = Effect + z_90 * se,
    cil_90 = Effect - z_90 * se
  )

cit <- c("Whataboutism\nvs. Identity Threat", lm2a$coefficients[2])
cit[2] <- lapply(cit[2], as.numeric); cit <- as.data.frame(cit)
colnames(cit) <- c("Treatment_Type", "Effect")
cit$se <- lm2a$std.error[2]
cit <- cit %>% 
  mutate(
    cih_95 = Effect + z_95 * se,
    cil_95 = Effect - z_95 * se,
    cih_90 = Effect + z_90 * se,
    cil_90 = Effect - z_90 * se
  )

com <- c("Whataboutism\nvs. Control", lm3a$coefficients[2])
com[2] <- lapply(com[2], as.numeric); com <- as.data.frame(com)
colnames(com) <- c("Treatment_Type", "Effect")
com$se <- lm3a$std.error[2]
com <- com %>% 
  mutate(
    cih_95 = Effect + z_95 * se,
    cil_95 = Effect - z_95 * se,
    cih_90 = Effect + z_90 * se,
    cil_90 = Effect - z_90 * se
  )

#bind results by treatment type
results_left <- rbind(unlist(idt), unlist(cit), unlist(com))
colnames(results_left) <- c("Treatment_Type", "Effect", "se", "cih_95", "cil_95", "cih_90", "cil_90")
results_left <- data.frame(results_left)
results_left$Effect <- as.numeric(results_left$Effect)
results_left$cil_95 <- as.numeric(results_left$cil_95)
results_left$cih_95 <- as.numeric(results_left$cih_95)
results_left$cil_90 <- as.numeric(results_left$cil_90)
results_left$cih_90 <- as.numeric(results_left$cih_90)
results_left$HTE <- "Left"

#Right-wing
data1_right <- subset(data1, data1$lr_1 >= 6)
data2_right <- subset(data2, data2$lr_1 >= 6)
data3_right <- subset(data3, data3$lr_1 >= 6)

#models by treatment type
lm1a <- lm_robust(prej_direct_1 ~ as.factor(trt) + education + income + newjob, data=data1_right, se_type="HC3")
lm2a <- lm_robust(prej_direct_1 ~ as.factor(trt) + lr_1 + natid_covar1_1 + income + newjob, data = data2_right, se_type="HC3")
lm3a <- lm_robust(prej_direct_1 ~ as.factor(trt) + lr_1 + unemployed, data = data3_right, se_type="HC3")

idt <- c("Identity Threat\nvs. Control", lm1a$coefficients[2])
idt[2] <- lapply(idt[2], as.numeric); idt <- as.data.frame(idt)
colnames(idt) <- c("Treatment_Type", "Effect")
idt$se <- lm1a$std.error[2]
idt <- idt %>% 
  mutate(
    cih_95 = Effect + z_95 * se,
    cil_95 = Effect - z_95 * se,
    cih_90 = Effect + z_90 * se,
    cil_90 = Effect - z_90 * se
  )

cit <- c("Whataboutism\nvs. Identity Threat", lm2a$coefficients[2])
cit[2] <- lapply(cit[2], as.numeric); cit <- as.data.frame(cit)
colnames(cit) <- c("Treatment_Type", "Effect")
cit$se <- lm2a$std.error[2]
cit <- cit %>% 
  mutate(
    cih_95 = Effect + z_95 * se,
    cil_95 = Effect - z_95 * se,
    cih_90 = Effect + z_90 * se,
    cil_90 = Effect - z_90 * se
  )

com <- c("Whataboutism\nvs. Control", lm3a$coefficients[2])
com[2] <- lapply(com[2], as.numeric); com <- as.data.frame(com)
colnames(com) <- c("Treatment_Type", "Effect")
com$se <- lm3a$std.error[2]
com <- com %>% 
  mutate(
    cih_95 = Effect + z_95 * se,
    cil_95 = Effect - z_95 * se,
    cih_90 = Effect + z_90 * se,
    cil_90 = Effect - z_90 * se
  )

#bind results by treatment type
results_right <- rbind(unlist(idt), unlist(cit), unlist(com))
colnames(results_right) <- c("Treatment_Type", "Effect", "se", "cih_95", "cil_95", "cih_90", "cil_90")
results_right <- data.frame(results_right)
results_right$Effect <- as.numeric(results_right$Effect)
results_right$cil_95 <- as.numeric(results_right$cil_95)
results_right$cih_95 <- as.numeric(results_right$cih_95)
results_right$cil_90 <- as.numeric(results_right$cil_90)
results_right$cih_90 <- as.numeric(results_right$cih_90)
results_right$HTE <- "Right"

#bind left and right results and plot
lr_hte <- rbind(results_left, results_right)
lr_hte$HTE <- factor(lr_hte$HTE, levels = c("Left", "Right"))

ep1 <- ggplot(lr_hte, aes(x = HTE, y = Effect, colour=factor(Treatment_Type, level=factor_levels))) +
  geom_point(stat="identity", position=position_dodge(width=0.3), size=3) + 
  labs(tag = "A") +
  scale_y_continuous(limits = c(-1.3, 2.1)) +
  geom_linerange(aes(ymin = cil_95, ymax = cih_95), 
                 position = position_dodge(width=0.3,),
                 size = 1) +
  geom_linerange(aes(ymin = cil_90, ymax = cih_90), 
                 position = position_dodge(width=0.3),
                 size = 2)+
  ylab("Treatment Effect") +
  xlab("Left-Right Self-Placement") +
  geom_hline(yintercept=0, linetype="dashed", color="blue") +
  labs(color='Update Type') +
  theme_bw()+
  theme(
    axis.text=element_text(size=16),
    plot.title = element_text(size=25),
    axis.title=element_text(size=18),
    legend.text=element_text(size=12),
    legend.title=element_text(size=16),
    legend.key.height = unit(1.5, 'cm')
  )

####TABLE D1(B) - Explicit Prejudice by Age####

#age group 1

data1_age1 <- subset(data1, data1$age<31)
data2_age1 <- subset(data2, data2$age<31)
data3_age1 <- subset(data3, data3$age<31)

#models by treatment type
lm1a <- lm_robust(prej_direct_1 ~ as.factor(trt) + education + income + newjob, data=data1_age1, se_type="HC3")
lm2a <- lm_robust(prej_direct_1 ~ as.factor(trt) + lr_1  + natid_covar1_1 + income + newjob, data=data2_age1, se_type="HC3")
lm3a <- lm_robust(prej_direct_1 ~ as.factor(trt) + lr_1 + unemployed, data=data3_age1, se_type="HC3")

idt <- c("Identity Threat\nvs. Control", lm1a$coefficients[2]) 
idt[2] <- lapply(idt[2], as.numeric); idt <- as.data.frame(idt)
colnames(idt) <- c("Treatment_Type", "Effect")
idt$se <- (lm1a$std.error)[2]
idt <- idt %>% mutate(cih90 = Effect + 1.645*se)
idt <- idt %>% mutate(cil90 = Effect - 1.645*se)
idt <- idt %>% mutate(cih95 = Effect + 1.96*se)
idt <- idt %>% mutate(cil95 = Effect - 1.96*se)

cit <- c("Whataboutism\nvs. Identity Threat", lm2a$coefficients[2]) 
cit[2] <- lapply(cit[2], as.numeric); cit <- as.data.frame(cit)
colnames(cit) <- c("Treatment_Type", "Effect")
cit$se <- (lm2a$std.error)[2]
cit <- cit %>% mutate(cih90 = Effect + 1.645*se)
cit <- cit %>% mutate(cil90 = Effect - 1.645*se)
cit <- cit %>% mutate(cih95 = Effect + 1.96*se)
cit <- cit %>% mutate(cil95 = Effect - 1.96*se)

com <- c("Whataboutism\nvs. Control", lm3a$coefficients[2]) 
com[2] <- lapply(com[2], as.numeric); com <- as.data.frame(com)
colnames(com) <- c("Treatment_Type", "Effect")
com$se <- (lm3a$std.error)[2]
com <- com %>% mutate(cih90 = Effect + 1.645*se)
com <- com %>% mutate(cil90 = Effect - 1.645*se)
com <- com %>% mutate(cih95 = Effect + 1.96*se)
com <- com %>% mutate(cil95 = Effect - 1.96*se)

#bind results by treatment type
results_age1 <- rbind(unlist(idt), unlist(cit), unlist(com))
colnames(results_age1) <- c("Treatment_Type", "Effect", "se", "cih90", "cil90", "cih95", "cil95")
results_age1 <- data.frame(results_age1)
results_age1$Effect <- as.numeric(results_age1$Effect)
results_age1$cil90 <- as.numeric(results_age1$cil90)
results_age1$cih90 <- as.numeric(results_age1$cih90)
results_age1$cil95 <- as.numeric(results_age1$cil95)
results_age1$cih95 <- as.numeric(results_age1$cih95)

results_age1$Age<- "18-30"

#age group 2

data1_age2 <- subset(data1, data1$age>30 & data1$age<41)
data2_age2 <- subset(data2, data2$age>30 & data2$age<41)
data3_age2 <- subset(data3, data3$age>30 & data3$age<41)

#models by treatment type
lm1a <- lm_robust(prej_direct_1 ~ as.factor(trt) + education + income + newjob, data=data1_age2, se_type="HC3")
lm2a <- lm_robust(prej_direct_1 ~ as.factor(trt) + lr_1  + natid_covar1_1 + income + newjob, data=data2_age2, se_type="HC3")
lm3a <- lm_robust(prej_direct_1 ~ as.factor(trt) + lr_1 + unemployed, data=data3_age2, se_type="HC3")

idt <- c("Identity Threat\nvs. Control", lm1a$coefficients[2]) 
idt[2] <- lapply(idt[2], as.numeric); idt <- as.data.frame(idt)
colnames(idt) <- c("Treatment_Type", "Effect")
idt$se <- (lm1a$std.error)[2]
idt <- idt %>% mutate(cih90 = Effect + 1.645*se)
idt <- idt %>% mutate(cil90 = Effect - 1.645*se)
idt <- idt %>% mutate(cih95 = Effect + 1.96*se)
idt <- idt %>% mutate(cil95 = Effect - 1.96*se)

cit <- c("Whataboutism\nvs. Identity Threat", lm2a$coefficients[2]) 
cit[2] <- lapply(cit[2], as.numeric); cit <- as.data.frame(cit)
colnames(cit) <- c("Treatment_Type", "Effect")
cit$se <- (lm2a$std.error)[2]
cit <- cit %>% mutate(cih90 = Effect + 1.645*se)
cit <- cit %>% mutate(cil90 = Effect - 1.645*se)
cit <- cit %>% mutate(cih95 = Effect + 1.96*se)
cit <- cit %>% mutate(cil95 = Effect - 1.96*se)

com <- c("Whataboutism\nvs. Control", lm3a$coefficients[2]) 
com[2] <- lapply(com[2], as.numeric); com <- as.data.frame(com)
colnames(com) <- c("Treatment_Type", "Effect")
com$se <- (lm3a$std.error)[2]
com <- com %>% mutate(cih90 = Effect + 1.645*se)
com <- com %>% mutate(cil90 = Effect - 1.645*se)
com <- com %>% mutate(cih95 = Effect + 1.96*se)
com <- com %>% mutate(cil95 = Effect - 1.96*se)

#bind results by treatment type
results_age2 <- rbind(unlist(idt), unlist(cit), unlist(com))
colnames(results_age2) <- c("Treatment_Type", "Effect", "se", "cih90", "cil90", "cih95", "cil95")
results_age2 <- data.frame(results_age2)
results_age2$Effect <- as.numeric(results_age2$Effect)
results_age2$cil90 <- as.numeric(results_age2$cil90)
results_age2$cih90 <- as.numeric(results_age2$cih90)
results_age2$cil95 <- as.numeric(results_age2$cil95)
results_age2$cih95 <- as.numeric(results_age2$cih95)
results_age2$Age<- "31-40"

#age group 3

data1_age3 <- subset(data1, data1$age>40 & data1$age<51)
data2_age3 <- subset(data2, data2$age>40 & data2$age<51)
data3_age3 <- subset(data3, data3$age>40 & data3$age<51)

#models by treatment type
lm1a <- lm_robust(prej_direct_1 ~ as.factor(trt) + education + income + newjob, data=data1_age3, se_type="HC3")
lm2a <- lm_robust(prej_direct_1 ~ as.factor(trt) + lr_1  + natid_covar1_1 + income + newjob, data=data2_age3, se_type="HC3")
lm3a <- lm_robust(prej_direct_1 ~ as.factor(trt) + lr_1 + unemployed, data=data3_age3, se_type="HC3")

idt <- c("Identity Threat\nvs. Control", lm1a$coefficients[2]) 
idt[2] <- lapply(idt[2], as.numeric); idt <- as.data.frame(idt)
colnames(idt) <- c("Treatment_Type", "Effect")
idt$se <- (lm1a$std.error)[2]
idt <- idt %>% mutate(cih90 = Effect + 1.645*se)
idt <- idt %>% mutate(cil90 = Effect - 1.645*se)
idt <- idt %>% mutate(cih95 = Effect + 1.96*se)
idt <- idt %>% mutate(cil95 = Effect - 1.96*se)

cit <- c("Whataboutism\nvs. Identity Threat", lm2a$coefficients[2]) 
cit[2] <- lapply(cit[2], as.numeric); cit <- as.data.frame(cit)
colnames(cit) <- c("Treatment_Type", "Effect")
cit$se <- (lm2a$std.error)[2]
cit <- cit %>% mutate(cih90 = Effect + 1.645*se)
cit <- cit %>% mutate(cil90 = Effect - 1.645*se)
cit <- cit %>% mutate(cih95 = Effect + 1.96*se)
cit <- cit %>% mutate(cil95 = Effect - 1.96*se)

com <- c("Whataboutism\nvs. Control", lm3a$coefficients[2]) 
com[2] <- lapply(com[2], as.numeric); com <- as.data.frame(com)
colnames(com) <- c("Treatment_Type", "Effect")
com$se <- (lm3a$std.error)[2]
com <- com %>% mutate(cih90 = Effect + 1.645*se)
com <- com %>% mutate(cil90 = Effect - 1.645*se)
com <- com %>% mutate(cih95 = Effect + 1.96*se)
com <- com %>% mutate(cil95 = Effect - 1.96*se)

#bind results by treatment type
results_age3 <- rbind(unlist(idt), unlist(cit), unlist(com))
colnames(results_age3) <- c("Treatment_Type", "Effect", "se", "cih90", "cil90", "cih95", "cil95")
results_age3 <- data.frame(results_age3)
results_age3$Effect <- as.numeric(results_age3$Effect)
results_age3$cil90 <- as.numeric(results_age3$cil90)
results_age3$cih90 <- as.numeric(results_age3$cih90)
results_age3$cil95 <- as.numeric(results_age3$cil95)
results_age3$cih95 <- as.numeric(results_age3$cih95)
results_age3$Age<- "41-50"

#age group 4

data1_age4 <- subset(data1, data1$age>50 & data1$age<61)
data2_age4 <- subset(data2, data2$age>50 & data2$age<61)
data3_age4 <- subset(data3, data3$age>50 & data3$age<61)

#models by treatment type
lm1a <- lm_robust(prej_direct_1 ~ as.factor(trt) + education + income + newjob, data=data1_age4, se_type="HC3")
lm2a <- lm_robust(prej_direct_1 ~ as.factor(trt) + lr_1  + natid_covar1_1 + income + newjob, data=data2_age4, se_type="HC3")
lm3a <- lm_robust(prej_direct_1 ~ as.factor(trt) + lr_1 + unemployed, data=data3_age4, se_type="HC3")

idt <- c("Identity Threat\nvs. Control", lm1a$coefficients[2]) 
idt[2] <- lapply(idt[2], as.numeric); idt <- as.data.frame(idt)
colnames(idt) <- c("Treatment_Type", "Effect")
idt$se <- (lm1a$std.error)[2]
idt <- idt %>% mutate(cih90 = Effect + 1.645*se)
idt <- idt %>% mutate(cil90 = Effect - 1.645*se)
idt <- idt %>% mutate(cih95 = Effect + 1.96*se)
idt <- idt %>% mutate(cil95 = Effect - 1.96*se)

cit <- c("Whataboutism\nvs. Identity Threat", lm2a$coefficients[2]) 
cit[2] <- lapply(cit[2], as.numeric); cit <- as.data.frame(cit)
colnames(cit) <- c("Treatment_Type", "Effect")
cit$se <- (lm2a$std.error)[2]
cit <- cit %>% mutate(cih90 = Effect + 1.645*se)
cit <- cit %>% mutate(cil90 = Effect - 1.645*se)
cit <- cit %>% mutate(cih95 = Effect + 1.96*se)
cit <- cit %>% mutate(cil95 = Effect - 1.96*se)

com <- c("Whataboutism\nvs. Control", lm3a$coefficients[2]) 
com[2] <- lapply(com[2], as.numeric); com <- as.data.frame(com)
colnames(com) <- c("Treatment_Type", "Effect")
com$se <- (lm3a$std.error)[2]
com <- com %>% mutate(cih90 = Effect + 1.645*se)
com <- com %>% mutate(cil90 = Effect - 1.645*se)
com <- com %>% mutate(cih95 = Effect + 1.96*se)
com <- com %>% mutate(cil95 = Effect - 1.96*se)

#bind results by treatment type
results_age4 <- rbind(unlist(idt), unlist(cit), unlist(com))
colnames(results_age4) <- c("Treatment_Type", "Effect", "se", "cih90", "cil90", "cih95", "cil95")
results_age4 <- data.frame(results_age4)
results_age4$Effect <- as.numeric(results_age4$Effect)
results_age4$cil90 <- as.numeric(results_age4$cil90)
results_age4$cih90 <- as.numeric(results_age4$cih90)
results_age4$cil95 <- as.numeric(results_age4$cil95)
results_age4$cih95 <- as.numeric(results_age4$cih95)
results_age4$Age<- "51-60"

#age group 5

data1_age5 <- subset(data1, data1$age>60)
data2_age5 <- subset(data2, data2$age>60)
data3_age5 <- subset(data3, data3$age>60)

#models by treatment type
lm1a <- lm_robust(prej_direct_1 ~ as.factor(trt) + education + income + newjob, data=data1_age5, se_type="HC3")
lm2a <- lm_robust(prej_direct_1 ~ as.factor(trt) + lr_1  + natid_covar1_1 + income + newjob, data=data2_age5, se_type="HC3")
lm3a <- lm_robust(prej_direct_1 ~ as.factor(trt) + lr_1 + unemployed, data=data3_age5, se_type="HC3")

idt <- c("Identity Threat\nvs. Control", lm1a$coefficients[2]) 
idt[2] <- lapply(idt[2], as.numeric); idt <- as.data.frame(idt)
colnames(idt) <- c("Treatment_Type", "Effect")
idt$se <- (lm1a$std.error)[2]
idt <- idt %>% mutate(cih90 = Effect + 1.645*se)
idt <- idt %>% mutate(cil90 = Effect - 1.645*se)
idt <- idt %>% mutate(cih95 = Effect + 1.96*se)
idt <- idt %>% mutate(cil95 = Effect - 1.96*se)

cit <- c("Whataboutism\nvs. Identity Threat", lm2a$coefficients[2]) 
cit[2] <- lapply(cit[2], as.numeric); cit <- as.data.frame(cit)
colnames(cit) <- c("Treatment_Type", "Effect")
cit$se <- (lm2a$std.error)[2]
cit <- cit %>% mutate(cih90 = Effect + 1.645*se)
cit <- cit %>% mutate(cil90 = Effect - 1.645*se)
cit <- cit %>% mutate(cih95 = Effect + 1.96*se)
cit <- cit %>% mutate(cil95 = Effect - 1.96*se)

com <- c("Whataboutism\nvs. Control", lm3a$coefficients[2]) 
com[2] <- lapply(com[2], as.numeric); com <- as.data.frame(com)
colnames(com) <- c("Treatment_Type", "Effect")
com$se <- (lm3a$std.error)[2]
com <- com %>% mutate(cih90 = Effect + 1.645*se)
com <- com %>% mutate(cil90 = Effect - 1.645*se)
com <- com %>% mutate(cih95 = Effect + 1.96*se)
com <- com %>% mutate(cil95 = Effect - 1.96*se)

#bind results by treatment type
results_age5 <- rbind(unlist(idt), unlist(cit), unlist(com))
colnames(results_age5) <- c("Treatment_Type", "Effect", "se", "cih90", "cil90", "cih95", "cil95")
results_age5 <- data.frame(results_age5)
results_age5$Effect <- as.numeric(results_age5$Effect)
results_age5$cil90 <- as.numeric(results_age5$cil90)
results_age5$cih90 <- as.numeric(results_age5$cih90)
results_age5$cil95 <- as.numeric(results_age5$cil95)
results_age5$cih95 <- as.numeric(results_age5$cih95)
results_age5$Age<- "61+"

#combine age groups and plot

age_hte <- rbind(results_age1, results_age2, results_age3, results_age4, results_age5)
age_hte$Age <- factor(age_hte$Age, levels = c("18-30", "31-40", "41-50", "51-60", "61+"))

ep2 <- ggplot(age_hte, aes(x = Age, y = Effect, colour=factor(Treatment_Type, level=factor_levels))) +
  geom_point(stat="identity", position=position_dodge(width=0.5), size=3) + 
  labs(tag = "B") +
  scale_y_continuous(limits = c(-1.3, 2.1)) +
  geom_linerange(aes(ymin = cil95, ymax = cih95), 
                 position = position_dodge(width=0.5,),
                 size = 1) +
  geom_linerange(aes(ymin = cil90, ymax = cih90), 
                 position = position_dodge(width=0.5),
                 size = 2) +
  ylab("Treatment Effect") +
  xlab("Age") +
  geom_hline(yintercept=0, linetype="dashed", color="blue") +
  labs(color='Update Type') +
  theme_bw()+
  theme(
    axis.text=element_text(size=16),
    plot.title = element_text(size=25),
    axis.title=element_text(size=18),
    legend.text=element_text(size=12),
    legend.title=element_text(size=16),
    legend.key.height = unit(1.5, 'cm')
  )

####TABLE D1(C) - Explicit Prejudice by Gender####

# women
data1_f <- subset(data1, data1$female == 1)
data2_f <- subset(data2, data2$female == 1)
data3_f <- subset(data3, data3$female == 1)

#models by treatment type
lm1a <- lm_robust(prej_direct_1 ~ as.factor(trt) + education + income + newjob, data=data1_f, se_type="HC3")
lm2a <- lm_robust(prej_direct_1 ~ as.factor(trt) + lr_1 + natid_covar1_1 + income + newjob, data = data2_f, se_type="HC3")
lm3a <- lm_robust(prej_direct_1 ~ as.factor(trt) + lr_1 + unemployed, data = data3_f, se_type="HC3")

idt <- data.frame(
  Treatment_Type = "Identity Threat\nvs. Control",
  Effect = lm1a$coefficients[2],
  se = lm1a$std.error[2]
)
idt <- idt %>% 
  mutate(
    ci90h = Effect + 1.645 * se,
    ci90l = Effect - 1.645 * se,
    ci95h = Effect + 1.96 * se,
    ci95l = Effect - 1.96 * se
  )

cit <- data.frame(
  Treatment_Type = "Whataboutism\nvs. Identity Threat",
  Effect = lm2a$coefficients[2],
  se = lm2a$std.error[2]
)
cit <- cit %>% 
  mutate(
    ci90h = Effect + 1.645 * se,
    ci90l = Effect - 1.645 * se,
    ci95h = Effect + 1.96 * se,
    ci95l = Effect - 1.96 * se
  )

com <- data.frame(
  Treatment_Type = "Whataboutism\nvs. Control",
  Effect = lm3a$coefficients[2],
  se = lm3a$std.error[2]
)
com <- com %>% 
  mutate(
    ci90h = Effect + 1.645 * se,
    ci90l = Effect - 1.645 * se,
    ci95h = Effect + 1.96 * se,
    ci95l = Effect - 1.96 * se
  )

#bind results by treatment type
results_f <- rbind(idt, cit, com)
results_f$Gender <- "Female"

# men
data1_m <- subset(data1, data1$female == 0)
data2_m <- subset(data2, data2$female == 0)
data3_m <- subset(data3, data3$female == 0)

#models by treatment type
lm1a <- lm_robust(prej_direct_1 ~ as.factor(trt) + education + income + newjob, data=data1_m, se_type="HC3")
lm2a <- lm_robust(prej_direct_1 ~ as.factor(trt) + lr_1 + natid_covar1_1 + income + newjob, data = data2_m, se_type="HC3")
lm3a <- lm_robust(prej_direct_1 ~ as.factor(trt) + lr_1 + unemployed, data = data3_m, se_type="HC3")

idt <- data.frame(
  Treatment_Type = "Identity Threat\nvs. Control",
  Effect = lm1a$coefficients[2],
  se = lm1a$std.error[2]
)
idt <- idt %>% 
  mutate(
    ci90h = Effect + 1.645 * se,
    ci90l = Effect - 1.645 * se,
    ci95h = Effect + 1.96 * se,
    ci95l = Effect - 1.96 * se
  )

cit <- data.frame(
  Treatment_Type = "Whataboutism\nvs. Identity Threat",
  Effect = lm2a$coefficients[2],
  se = lm2a$std.error[2]
)
cit <- cit %>% 
  mutate(
    ci90h = Effect + 1.645 * se,
    ci90l = Effect - 1.645 * se,
    ci95h = Effect + 1.96 * se,
    ci95l = Effect - 1.96 * se
  )

com <- data.frame(
  Treatment_Type = "Whataboutism\nvs. Control",
  Effect = lm3a$coefficients[2],
  se = lm3a$std.error[2]
)
com <- com %>% 
  mutate(
    ci90h = Effect + 1.645 * se,
    ci90l = Effect - 1.645 * se,
    ci95h = Effect + 1.96 * se,
    ci95l = Effect - 1.96 * se
  )

#bind results by treatment type
results_m <- rbind(idt, cit, com)
results_m$Gender <- "Male"

# Combine women/men and plot
gender_hte <- rbind(results_f, results_m)
gender_hte$Gender <- factor(gender_hte$Gender, levels = c("Female", "Male"))

ep3 <- ggplot(gender_hte, aes(x = Gender, y = Effect, colour=factor(Treatment_Type, level=factor_levels))) +
  geom_point(stat="identity", position=position_dodge(width=0.3), size=3) + 
  labs(tag = "C") +
  scale_y_continuous(limits = c(-1.3, 2.1)) +
  geom_linerange(aes(ymin = ci95l, ymax = ci95h), 
                 position = position_dodge(width=0.3,),
                 size = 1) +
  geom_linerange(aes(ymin = ci90l, ymax = ci90h), 
                 position = position_dodge(width=0.3),
                 size = 2)+
  ylab("Treatment Effect") +
  xlab("Gender") +
  geom_hline(yintercept=0, linetype="dashed", color="blue") +
  labs(color='Update Type') +
  theme_bw()+
  theme(
    axis.text=element_text(size=16),
    plot.title = element_text(size=25),
    axis.title=element_text(size=18),
    legend.text=element_text(size=12),
    legend.title=element_text(size=16),
    legend.key.height = unit(1.5, 'cm')
  )


####TABLE D1(D) - Explicit Prejudice by Pre-Treatment National Identification####

#quintile 1

data1_nid1 <- subset(data1, data1$nat_quint==1)
data2_nid1 <- subset(data2, data2$nat_quint==1)
data3_nid1 <- subset(data3, data3$nat_quint==1)

#models by treatment type
lm1a <- lm_robust(prej_direct_1 ~ as.factor(trt) + education + income + newjob, data=data1_nid1, se_type="HC3")
lm2a <- lm_robust(prej_direct_1 ~ as.factor(trt) + lr_1 + natid_covar1_1 + income + newjob, data = data2_nid1, se_type="HC3")
lm3a <- lm_robust(prej_direct_1 ~ as.factor(trt) + lr_1 + unemployed, data = data3_nid1, se_type="HC3")

idt <- c("Identity Threat\nvs. Control", lm1a$coefficients[2]) 
idt[2] <- lapply(idt[2], as.numeric); idt <- as.data.frame(idt)
colnames(idt) <- c("Treatment_Type", "Effect")
idt$se <- (lm1a$std.error)[2]
idt <- idt %>% mutate(cih95 = Effect + 1.96*se)
idt <- idt %>% mutate(cil95 = Effect - 1.96*se)
idt <- idt %>% mutate(cih90 = Effect + 1.645*se)
idt <- idt %>% mutate(cil90 = Effect - 1.645*se)

cit <- c("Whataboutism\nvs. Identity Threat", lm2a$coefficients[2]) 
cit[2] <- lapply(cit[2], as.numeric); cit <- as.data.frame(cit)
colnames(cit) <- c("Treatment_Type", "Effect")
cit$se <- (lm2a$std.error)[2]
cit <- cit %>% mutate(cih95 = Effect + 1.96*se)
cit <- cit %>% mutate(cil95 = Effect - 1.96*se)
cit <- cit %>% mutate(cih90 = Effect + 1.645*se)
cit <- cit %>% mutate(cil90 = Effect - 1.645*se)

com <- c("Whataboutism\nvs. Control", lm3a$coefficients[2]) 
com[2] <- lapply(com[2], as.numeric); com <- as.data.frame(com)
colnames(com) <- c("Treatment_Type", "Effect")
com$se <- (lm3a$std.error)[2]
com <- com %>% mutate(cih95 = Effect + 1.96*se)
com <- com %>% mutate(cil95 = Effect - 1.96*se)
com <- com %>% mutate(cih90 = Effect + 1.645*se)
com <- com %>% mutate(cil90 = Effect - 1.645*se)

#bind results by treatment type
results_nid1 <- rbind(unlist(idt), unlist(cit), unlist(com))
colnames(results_nid1) <- c("Treatment_Type", "Effect", "se", "cih95", "cil95", "cih90", "cil90")
results_nid1 <- data.frame(results_nid1)
results_nid1$Effect <- as.numeric(results_nid1$Effect)
results_nid1$cil95 <- as.numeric(results_nid1$cil95)
results_nid1$cil90 <- as.numeric(results_nid1$cil90)
results_nid1$cih95 <- as.numeric(results_nid1$cih95)
results_nid1$cih90 <- as.numeric(results_nid1$cih90)
results_nid1$Quintile<- 1

#quintile 2

data1_nid2 <- subset(data1, data1$nat_quint==2)
data2_nid2 <- subset(data2, data2$nat_quint==2)
data3_nid2 <- subset(data3, data3$nat_quint==2)

#models by treatment type
lm1a <- lm_robust(prej_direct_1 ~ as.factor(trt) + education + income + newjob, data=data1_nid2, se_type="HC3")
lm2a <- lm_robust(prej_direct_1 ~ as.factor(trt) + lr_1 + natid_covar1_1 + income + newjob, data = data2_nid2, se_type="HC3")
lm3a <- lm_robust(prej_direct_1 ~ as.factor(trt) + lr_1 + unemployed, data = data3_nid2, se_type="HC3")

idt <- c("Identity Threat\nvs. Control", lm1a$coefficients[2]) 
idt[2] <- lapply(idt[2], as.numeric); idt <- as.data.frame(idt)
colnames(idt) <- c("Treatment_Type", "Effect")
idt$se <- (lm1a$std.error)[2]
idt <- idt %>% mutate(cih95 = Effect + 1.96*se)
idt <- idt %>% mutate(cil95 = Effect - 1.96*se)
idt <- idt %>% mutate(cih90 = Effect + 1.645*se)
idt <- idt %>% mutate(cil90 = Effect - 1.645*se)

cit <- c("Whataboutism\nvs. Identity Threat", lm2a$coefficients[2]) 
cit[2] <- lapply(cit[2], as.numeric); cit <- as.data.frame(cit)
colnames(cit) <- c("Treatment_Type", "Effect")
cit$se <- (lm2a$std.error)[2]
cit <- cit %>% mutate(cih95 = Effect + 1.96*se)
cit <- cit %>% mutate(cil95 = Effect - 1.96*se)
cit <- cit %>% mutate(cih90 = Effect + 1.645*se)
cit <- cit %>% mutate(cil90 = Effect - 1.645*se)

com <- c("Whataboutism\nvs. Control", lm3a$coefficients[2]) 
com[2] <- lapply(com[2], as.numeric); com <- as.data.frame(com)
colnames(com) <- c("Treatment_Type", "Effect")
com$se <- (lm3a$std.error)[2]
com <- com %>% mutate(cih95 = Effect + 1.96*se)
com <- com %>% mutate(cil95 = Effect - 1.96*se)
com <- com %>% mutate(cih90 = Effect + 1.645*se)
com <- com %>% mutate(cil90 = Effect - 1.645*se)

#bind results by treatment type
results_nid2 <- rbind(unlist(idt), unlist(cit), unlist(com))
colnames(results_nid2) <- c("Treatment_Type", "Effect", "se", "cih95", "cil95", "cih90", "cil90")
results_nid2 <- data.frame(results_nid2)
results_nid2$Effect <- as.numeric(results_nid2$Effect)
results_nid2$cil95 <- as.numeric(results_nid2$cil95)
results_nid2$cil90 <- as.numeric(results_nid2$cil90)
results_nid2$cih95 <- as.numeric(results_nid2$cih95)
results_nid2$cih90 <- as.numeric(results_nid2$cih90)
results_nid2$Quintile<- 2

#quintile 3

data1_nid3 <- subset(data1, data1$nat_quint==3)
data2_nid3 <- subset(data2, data2$nat_quint==3)
data3_nid3 <- subset(data3, data3$nat_quint==3)

#models by treatment type
lm1a <- lm_robust(prej_direct_1 ~ as.factor(trt) + education + income + newjob, data=data1_nid3, se_type="HC3")
lm2a <- lm_robust(prej_direct_1 ~ as.factor(trt) + lr_1 + natid_covar1_1 + income + newjob, data = data2_nid3, se_type="HC3")
lm3a <- lm_robust(prej_direct_1 ~ as.factor(trt) + lr_1 + unemployed, data = data3_nid3, se_type="HC3")

idt <- c("Identity Threat\nvs. Control", lm1a$coefficients[2]) 
idt[2] <- lapply(idt[2], as.numeric); idt <- as.data.frame(idt)
colnames(idt) <- c("Treatment_Type", "Effect")
idt$se <- (lm1a$std.error)[2]
idt <- idt %>% mutate(cih95 = Effect + 1.96*se)
idt <- idt %>% mutate(cil95 = Effect - 1.96*se)
idt <- idt %>% mutate(cih90 = Effect + 1.645*se)
idt <- idt %>% mutate(cil90 = Effect - 1.645*se)

cit <- c("Whataboutism\nvs. Identity Threat", lm2a$coefficients[2]) 
cit[2] <- lapply(cit[2], as.numeric); cit <- as.data.frame(cit)
colnames(cit) <- c("Treatment_Type", "Effect")
cit$se <- (lm2a$std.error)[2]
cit <- cit %>% mutate(cih95 = Effect + 1.96*se)
cit <- cit %>% mutate(cil95 = Effect - 1.96*se)
cit <- cit %>% mutate(cih90 = Effect + 1.645*se)
cit <- cit %>% mutate(cil90 = Effect - 1.645*se)

com <- c("Whataboutism\nvs. Control", lm3a$coefficients[2]) 
com[2] <- lapply(com[2], as.numeric); com <- as.data.frame(com)
colnames(com) <- c("Treatment_Type", "Effect")
com$se <- (lm3a$std.error)[2]
com <- com %>% mutate(cih95 = Effect + 1.96*se)
com <- com %>% mutate(cil95 = Effect - 1.96*se)
com <- com %>% mutate(cih90 = Effect + 1.645*se)
com <- com %>% mutate(cil90 = Effect - 1.645*se)

#bind results by treatment type
results_nid3 <- rbind(unlist(idt), unlist(cit), unlist(com))
colnames(results_nid3) <- c("Treatment_Type", "Effect", "se", "cih95", "cil95", "cih90", "cil90")
results_nid3 <- data.frame(results_nid3)
results_nid3$Effect <- as.numeric(results_nid3$Effect)
results_nid3$cil95 <- as.numeric(results_nid3$cil95)
results_nid3$cil90 <- as.numeric(results_nid3$cil90)
results_nid3$cih95 <- as.numeric(results_nid3$cih95)
results_nid3$cih90 <- as.numeric(results_nid3$cih90)
results_nid3$Quintile<- 3

#quintile 4

data1_nid4 <- subset(data1, data1$nat_quint==4)
data2_nid4 <- subset(data2, data2$nat_quint==4)
data3_nid4 <- subset(data3, data3$nat_quint==4)

#models by treatment type
lm1a <- lm_robust(prej_direct_1 ~ as.factor(trt) + education + income + newjob, data=data1_nid4, se_type="HC3")
lm2a <- lm_robust(prej_direct_1 ~ as.factor(trt) + lr_1 + natid_covar1_1 + income + newjob, data = data2_nid4, se_type="HC3")
lm3a <- lm_robust(prej_direct_1 ~ as.factor(trt) + lr_1 + unemployed, data = data3_nid4, se_type="HC3")

idt <- c("Identity Threat\nvs. Control", lm1a$coefficients[2]) 
idt[2] <- lapply(idt[2], as.numeric); idt <- as.data.frame(idt)
colnames(idt) <- c("Treatment_Type", "Effect")
idt$se <- (lm1a$std.error)[2]
idt <- idt %>% mutate(cih95 = Effect + 1.96*se)
idt <- idt %>% mutate(cil95 = Effect - 1.96*se)
idt <- idt %>% mutate(cih90 = Effect + 1.645*se)
idt <- idt %>% mutate(cil90 = Effect - 1.645*se)

cit <- c("Whataboutism\nvs. Identity Threat", lm2a$coefficients[2]) 
cit[2] <- lapply(cit[2], as.numeric); cit <- as.data.frame(cit)
colnames(cit) <- c("Treatment_Type", "Effect")
cit$se <- (lm2a$std.error)[2]
cit <- cit %>% mutate(cih95 = Effect + 1.96*se)
cit <- cit %>% mutate(cil95 = Effect - 1.96*se)
cit <- cit %>% mutate(cih90 = Effect + 1.645*se)
cit <- cit %>% mutate(cil90 = Effect - 1.645*se)

com <- c("Whataboutism\nvs. Control", lm3a$coefficients[2]) 
com[2] <- lapply(com[2], as.numeric); com <- as.data.frame(com)
colnames(com) <- c("Treatment_Type", "Effect")
com$se <- (lm3a$std.error)[2]
com <- com %>% mutate(cih95 = Effect + 1.96*se)
com <- com %>% mutate(cil95 = Effect - 1.96*se)
com <- com %>% mutate(cih90 = Effect + 1.645*se)
com <- com %>% mutate(cil90 = Effect - 1.645*se)

#bind results by treatment type
results_nid4 <- rbind(unlist(idt), unlist(cit), unlist(com))
colnames(results_nid4) <- c("Treatment_Type", "Effect", "se", "cih95", "cil95", "cih90", "cil90")
results_nid4 <- data.frame(results_nid4)
results_nid4$Effect <- as.numeric(results_nid4$Effect)
results_nid4$cil95 <- as.numeric(results_nid4$cil95)
results_nid4$cil90 <- as.numeric(results_nid4$cil90)
results_nid4$cih95 <- as.numeric(results_nid4$cih95)
results_nid4$cih90 <- as.numeric(results_nid4$cih90)
results_nid4$Quintile<- 4

#quintile 5

data1_nid5 <- subset(data1, data1$nat_quint==5)
data2_nid5 <- subset(data2, data2$nat_quint==5)
data3_nid5 <- subset(data3, data3$nat_quint==5)

#models by treatment type
lm1a <- lm_robust(prej_direct_1 ~ as.factor(trt) + education + income + newjob, data=data1_nid5, se_type="HC3")
lm2a <- lm_robust(prej_direct_1 ~ as.factor(trt) + lr_1 + natid_covar1_1 + income + newjob, data = data2_nid5, se_type="HC3")
lm3a <- lm_robust(prej_direct_1 ~ as.factor(trt) + lr_1 + unemployed, data = data3_nid5, se_type="HC3")

idt <- c("Identity Threat\nvs. Control", lm1a$coefficients[2]) 
idt[2] <- lapply(idt[2], as.numeric); idt <- as.data.frame(idt)
colnames(idt) <- c("Treatment_Type", "Effect")
idt$se <- (lm1a$std.error)[2]
idt <- idt %>% mutate(cih95 = Effect + 1.96*se)
idt <- idt %>% mutate(cil95 = Effect - 1.96*se)
idt <- idt %>% mutate(cih90 = Effect + 1.645*se)
idt <- idt %>% mutate(cil90 = Effect - 1.645*se)

cit <- c("Whataboutism\nvs. Identity Threat", lm2a$coefficients[2]) 
cit[2] <- lapply(cit[2], as.numeric); cit <- as.data.frame(cit)
colnames(cit) <- c("Treatment_Type", "Effect")
cit$se <- (lm2a$std.error)[2]
cit <- cit %>% mutate(cih95 = Effect + 1.96*se)
cit <- cit %>% mutate(cil95 = Effect - 1.96*se)
cit <- cit %>% mutate(cih90 = Effect + 1.645*se)
cit <- cit %>% mutate(cil90 = Effect - 1.645*se)

com <- c("Whataboutism\nvs. Control", lm3a$coefficients[2]) 
com[2] <- lapply(com[2], as.numeric); com <- as.data.frame(com)
colnames(com) <- c("Treatment_Type", "Effect")
com$se <- (lm3a$std.error)[2]
com <- com %>% mutate(cih95 = Effect + 1.96*se)
com <- com %>% mutate(cil95 = Effect - 1.96*se)
com <- com %>% mutate(cih90 = Effect + 1.645*se)
com <- com %>% mutate(cil90 = Effect - 1.645*se)

#bind results by treatment type
results_nid5 <- rbind(unlist(idt), unlist(cit), unlist(com))
colnames(results_nid5) <- c("Treatment_Type", "Effect", "se", "cih95", "cil95", "cih90", "cil90")
results_nid5 <- data.frame(results_nid5)
results_nid5$Effect <- as.numeric(results_nid5$Effect)
results_nid5$cil95 <- as.numeric(results_nid5$cil95)
results_nid5$cil90 <- as.numeric(results_nid5$cil90)
results_nid5$cih95 <- as.numeric(results_nid5$cih95)
results_nid5$cih90 <- as.numeric(results_nid5$cih90)
results_nid5$Quintile<- 5

#combine quintiles and plot
nid_hte <- rbind(results_nid1, results_nid2, results_nid3, results_nid4, results_nid5)
nid_hte$Effect <- as.numeric(as.character(nid_hte$Effect))

ep4 <- ggplot(nid_hte, aes(x = Quintile, y = Effect, colour=factor(Treatment_Type, level=factor_levels))) +
  geom_point(stat="identity", position=position_dodge(width=0.5), size=3) + 
  labs(tag = "D") +
  scale_y_continuous(limits = c(-1.3, 2.1)) +
  geom_linerange(aes(ymin = cil95, ymax = cih95), 
                 position = position_dodge(width=0.5,),
                 size = 1) +
  geom_linerange(aes(ymin = cil90, ymax = cih90), 
                 position = position_dodge(width=0.5),
                 size = 2)+
  ylab("Treatment Effect") +
  xlab("National Identification") +
  geom_hline(yintercept=0, linetype="dashed", color="blue") +
  labs(color='Update Type') +
  theme_bw()+
  theme(
    axis.text=element_text(size=16),
    plot.title = element_text(size=25),
    axis.title=element_text(size=18),
    legend.text=element_text(size=12),
    legend.title=element_text(size=16),
    legend.key.height = unit(1.5, 'cm')
  )

####Combined Panel for Figure D1####

(ep1 | ep2) /
  (ep3 | ep4)+
  plot_layout(guides = "collect") &
  theme(legend.position = "right")

ggsave("figD1.png", width = 12, height = 10, dpi = 300)

####TABLE D2 - HTEs for Implicit Prejudice####

####TABLE D2(A) - Implicit Prejudice by Left/Right Self-Placement####

#left-wing

#models by treatment type
lm1i <- lm_robust(list_items ~ as.factor(trt) + education + income + newjob + as.factor(trt) * sens, data=data1_left, se_type="HC3")
lm2i <- lm_robust(list_items ~ as.factor(trt) + lr_1 + natid_covar1_1 + income + newjob + as.factor(trt) * sens, data=data2_left, se_type="HC3")
lm3i <- lm_robust(list_items ~ as.factor(trt) + lr_1 + unemployed + as.factor(trt) * sens, data=data3_left, se_type="HC3")

# calculate confidence intervals
com1 <- c(
  "Identity Threat\nvs. Control", 
  lm1i$coefficients[7], 
  lm1i$coefficients[7] + z_95 * lm1i$std.error[7], 
  lm1i$coefficients[7] - z_95 * lm1i$std.error[7], 
  lm1i$coefficients[7] + z_90 * lm1i$std.error[7], 
  lm1i$coefficients[7] - z_90 * lm1i$std.error[7]
)
com2 <- c(
  "Whataboutism\nvs. Identity Threat", 
  lm2i$coefficients[8], 
  lm2i$coefficients[8] + z_95 * lm2i$std.error[8], 
  lm2i$coefficients[8] - z_95 * lm2i$std.error[8], 
  lm2i$coefficients[8] + z_90 * lm2i$std.error[8], 
  lm2i$coefficients[8] - z_90 * lm2i$std.error[8]
)
com3 <- c(
  "Whataboutism\nvs. Control", 
  lm3i$coefficients[6], 
  lm3i$coefficients[6] + z_95 * lm3i$std.error[6], 
  lm3i$coefficients[6] - z_95 * lm3i$std.error[6], 
  lm3i$coefficients[6] + z_90 * lm3i$std.error[6], 
  lm3i$coefficients[6] - z_90 * lm3i$std.error[6]
)

#bind results by treatment type
com_l <- data.frame(rbind(com1, com2, com3))
com_l$Self_placement <- "Left"

#right-wing

#models by treatment type
lm1i <- lm_robust(list_items ~ as.factor(trt) + education + income + newjob + as.factor(trt) * sens, data=data1_right, se_type="HC3")
lm2i <- lm_robust(list_items ~ as.factor(trt) + lr_1 + natid_covar1_1 + income + newjob + as.factor(trt) * sens, data=data2_right, se_type="HC3")
lm3i <- lm_robust(list_items ~ as.factor(trt) + lr_1 + unemployed + as.factor(trt) * sens, data=data3_right, se_type="HC3")

# calculate confidence intervals
com1 <- c(
  "Identity Threat\nvs. Control", 
  lm1i$coefficients[7], 
  lm1i$coefficients[7] + z_95 * lm1i$std.error[7], 
  lm1i$coefficients[7] - z_95 * lm1i$std.error[7], 
  lm1i$coefficients[7] + z_90 * lm1i$std.error[7], 
  lm1i$coefficients[7] - z_90 * lm1i$std.error[7]
)
com2 <- c(
  "Whataboutism\nvs. Identity Threat", 
  lm2i$coefficients[8], 
  lm2i$coefficients[8] + z_95 * lm2i$std.error[8], 
  lm2i$coefficients[8] - z_95 * lm2i$std.error[8], 
  lm2i$coefficients[8] + z_90 * lm2i$std.error[8], 
  lm2i$coefficients[8] - z_90 * lm2i$std.error[8]
)
com3 <- c(
  "Whataboutism\nvs. Control", 
  lm3i$coefficients[6], 
  lm3i$coefficients[6] + z_95 * lm3i$std.error[6], 
  lm3i$coefficients[6] - z_95 * lm3i$std.error[6], 
  lm3i$coefficients[6] + z_90 * lm3i$std.error[6], 
  lm3i$coefficients[6] - z_90 * lm3i$std.error[6]
)

#bind results by treatment type
com_r <- data.frame(rbind(com1, com2, com3))
com_r$Self_placement <- "Right"

# Combine left and right
lr_hte <- rbind(com_l, com_r)
colnames(lr_hte) <- c("treatment", "effect", "cih95", "cil95", "cih90", "cil90", "Self_Placement")

# Ensure numeric data types for relevant columns
lr_hte$effect <- as.numeric(lr_hte$effect)
lr_hte$cih95 <- as.numeric(lr_hte$cih95)
lr_hte$cil95 <- as.numeric(lr_hte$cil95)
lr_hte$cih90 <- as.numeric(lr_hte$cih90)
lr_hte$cil90 <- as.numeric(lr_hte$cil90)

#plot
ip1 <- ggplot(lr_hte, aes(x = Self_Placement, y = effect, colour=factor(treatment, level=factor_levels))) +
  geom_point(stat="identity", position=position_dodge(width=0.3), size=3) + 
  labs(tag = "A") +
  scale_y_continuous(limits = c(-0.8, 0.8)) +
  geom_linerange(aes(ymin = cil95, ymax = cih95), 
                 position = position_dodge(width=0.3,),
                 size = 1) +
  geom_linerange(aes(ymin = cil90, ymax = cih90), 
                 position = position_dodge(width=0.3),
                 size = 2) +
  ylab("Treatment Effect") +
  xlab("Left-Right Self-Placement") +
  geom_hline(yintercept=0, linetype="dashed", color="blue") +
  labs(color='Update Type') +
  theme_bw()+
  theme(
    axis.text=element_text(size=16),
    plot.title = element_text(size=25),
    axis.title=element_text(size=18),
    legend.text=element_text(size=12),
    legend.title=element_text(size=16),
    legend.key.height = unit(1.5, 'cm')
  )


####TABLE D2(B) - Implicit Prejudice by Age####

#age group 1

#models by treatment type
lm1i <- lm_robust(list_items ~ as.factor(trt) + education + income + newjob + as.factor(trt) * sens, data=data1_age1, se_type="HC3")
lm2i <- lm_robust(list_items ~ as.factor(trt) + lr_1 + natid_covar1_1 + income + newjob + as.factor(trt) * sens, data=data2_age1, se_type="HC3")
lm3i <- lm_robust(list_items ~ as.factor(trt) + lr_1 + unemployed + as.factor(trt) * sens, data=data3_age1, se_type="HC3")

idt <- c("Identity Threat\nvs. Control", lm1i$coefficients[7]) 
idt[2] <- lapply(idt[2], as.numeric); idt <- as.data.frame(idt)
colnames(idt) <- c("Treatment_Type", "Effect")
idt$se <- (lm1i$std.error)[7]
idt <- idt %>% mutate(cih90 = Effect + 1.645*se)
idt <- idt %>% mutate(cil90 = Effect - 1.645*se)
idt <- idt %>% mutate(cih95 = Effect + 1.96*se)
idt <- idt %>% mutate(cil95 = Effect - 1.96*se)

cit <- c("Whataboutism\nvs. Identity Threat", lm2i$coefficients[8]) 
cit[2] <- lapply(cit[2], as.numeric); cit <- as.data.frame(cit)
colnames(cit) <- c("Treatment_Type", "Effect")
cit$se <- (lm2i$std.error)[8]
cit <- cit %>% mutate(cih90 = Effect + 1.645*se)
cit <- cit %>% mutate(cil90 = Effect - 1.645*se)
cit <- cit %>% mutate(cih95 = Effect + 1.96*se)
cit <- cit %>% mutate(cil95 = Effect - 1.96*se)

com <- c("Whataboutism\nvs. Control", lm3i$coefficients[6]) 
com[2] <- lapply(com[2], as.numeric); com <- as.data.frame(com)
colnames(com) <- c("Treatment Type", "Effect")
com$se <- (lm3i$std.error)[6]
com <- com %>% mutate(cih90 = Effect + 1.645*se)
com <- com %>% mutate(cil90 = Effect - 1.645*se)
com <- com %>% mutate(cih95 = Effect + 1.96*se)
com <- com %>% mutate(cil95 = Effect - 1.96*se)

#bind results by treatment type
results_age1 <- rbind(unlist(idt), unlist(cit), unlist(com))
colnames(results_age1) <- c("Treatment_Type", "Effect", "se", "cih90", "cil90", "cih95", "cil95")
results_age1 <- data.frame(results_age1)
results_age1$Effect <- as.numeric(results_age1$Effect)
results_age1$cil90 <- as.numeric(results_age1$cil90)
results_age1$cih90 <- as.numeric(results_age1$cih90)
results_age1$cil95 <- as.numeric(results_age1$cil95)
results_age1$cih95 <- as.numeric(results_age1$cih95)
results_age1$Age<- "18-30"

#age group 2

lm1i <- lm_robust(list_items ~ as.factor(trt) + education + income + newjob + as.factor(trt) * sens, data=data1_age2, se_type="HC3")
lm2i <- lm_robust(list_items ~ as.factor(trt) + lr_1 + natid_covar1_1 + income + newjob + as.factor(trt) * sens, data=data2_age2, se_type="HC3")
lm3i <- lm_robust(list_items ~ as.factor(trt) + lr_1 + unemployed + as.factor(trt) * sens, data=data3_age2, se_type="HC3")

idt <- c("Identity Threat\nvs. Control", lm1i$coefficients[7]) 
idt[2] <- lapply(idt[2], as.numeric); idt <- as.data.frame(idt)
colnames(idt) <- c("Treatment_Type", "Effect")
idt$se <- (lm1i$std.error)[7]
idt <- idt %>% mutate(cih90 = Effect + 1.645*se)
idt <- idt %>% mutate(cil90 = Effect - 1.645*se)
idt <- idt %>% mutate(cih95 = Effect + 1.96*se)
idt <- idt %>% mutate(cil95 = Effect - 1.96*se)

cit <- c("Whataboutism\nvs. Identity Threat", lm2i$coefficients[8]) 
cit[2] <- lapply(cit[2], as.numeric); cit <- as.data.frame(cit)
colnames(cit) <- c("Treatment_Type", "Effect")
cit$se <- (lm2i$std.error)[8]
cit <- cit %>% mutate(cih90 = Effect + 1.645*se)
cit <- cit %>% mutate(cil90 = Effect - 1.645*se)
cit <- cit %>% mutate(cih95 = Effect + 1.96*se)
cit <- cit %>% mutate(cil95 = Effect - 1.96*se)

com <- c("Whataboutism\nvs. Control", lm3i$coefficients[6]) 
com[2] <- lapply(com[2], as.numeric); com <- as.data.frame(com)
colnames(com) <- c("Treatment Type", "Effect")
com$se <- (lm3i$std.error)[6]
com <- com %>% mutate(cih90 = Effect + 1.645*se)
com <- com %>% mutate(cil90 = Effect - 1.645*se)
com <- com %>% mutate(cih95 = Effect + 1.96*se)
com <- com %>% mutate(cil95 = Effect - 1.96*se)

#bind results by treatment type
results_age2 <- rbind(unlist(idt), unlist(cit), unlist(com))
colnames(results_age2) <- c("Treatment_Type", "Effect", "se", "cih90", "cil90", "cih95", "cil95")
results_age2 <- data.frame(results_age2)
results_age2$Effect <- as.numeric(results_age2$Effect)
results_age2$cil90 <- as.numeric(results_age2$cil90)
results_age2$cih90 <- as.numeric(results_age2$cih90)
results_age2$cil95 <- as.numeric(results_age2$cil95)
results_age2$cih95 <- as.numeric(results_age2$cih95)
results_age2$Age<- "31-40"

#age group 3

lm1i <- lm_robust(list_items ~ as.factor(trt) + education + income + newjob + as.factor(trt) * sens, data=data1_age3, se_type="HC3")
lm2i <- lm_robust(list_items ~ as.factor(trt) + lr_1 + natid_covar1_1 + income + newjob + as.factor(trt) * sens, data=data2_age3, se_type="HC3")
lm3i <- lm_robust(list_items ~ as.factor(trt) + lr_1 + unemployed + as.factor(trt) * sens, data=data3_age3, se_type="HC3")

idt <- c("Identity Threat\nvs. Control", lm1i$coefficients[7]) 
idt[2] <- lapply(idt[2], as.numeric); idt <- as.data.frame(idt)
colnames(idt) <- c("Treatment_Type", "Effect")
idt$se <- (lm1i$std.error)[7]
idt <- idt %>% mutate(cih90 = Effect + 1.645*se)
idt <- idt %>% mutate(cil90 = Effect - 1.645*se)
idt <- idt %>% mutate(cih95 = Effect + 1.96*se)
idt <- idt %>% mutate(cil95 = Effect - 1.96*se)

cit <- c("Whataboutism\nvs. Identity Threat", lm2i$coefficients[8]) 
cit[2] <- lapply(cit[2], as.numeric); cit <- as.data.frame(cit)
colnames(cit) <- c("Treatment_Type", "Effect")
cit$se <- (lm2i$std.error)[8]
cit <- cit %>% mutate(cih90 = Effect + 1.645*se)
cit <- cit %>% mutate(cil90 = Effect - 1.645*se)
cit <- cit %>% mutate(cih95 = Effect + 1.96*se)
cit <- cit %>% mutate(cil95 = Effect - 1.96*se)

com <- c("Whataboutism\nvs. Control", lm3i$coefficients[6]) 
com[2] <- lapply(com[2], as.numeric); com <- as.data.frame(com)
colnames(com) <- c("Treatment Type", "Effect")
com$se <- (lm3i$std.error)[6]
com <- com %>% mutate(cih90 = Effect + 1.645*se)
com <- com %>% mutate(cil90 = Effect - 1.645*se)
com <- com %>% mutate(cih95 = Effect + 1.96*se)
com <- com %>% mutate(cil95 = Effect - 1.96*se)

#bind results by treatment type
results_age3 <- rbind(unlist(idt), unlist(cit), unlist(com))
colnames(results_age3) <- c("Treatment_Type", "Effect", "se", "cih90", "cil90", "cih95", "cil95")
results_age3 <- data.frame(results_age3)
results_age3$Effect <- as.numeric(results_age3$Effect)
results_age3$cil90 <- as.numeric(results_age3$cil90)
results_age3$cih90 <- as.numeric(results_age3$cih90)
results_age3$cil95 <- as.numeric(results_age3$cil95)
results_age3$cih95 <- as.numeric(results_age3$cih95)
results_age3$Age<- "41-50"

#age group 4

lm1i <- lm_robust(list_items ~ as.factor(trt) + education + income + newjob + as.factor(trt) * sens, data=data1_age4, se_type="HC3")
lm2i <- lm_robust(list_items ~ as.factor(trt) + lr_1 + natid_covar1_1 + income + newjob + as.factor(trt) * sens, data=data2_age4, se_type="HC3")
lm3i <- lm_robust(list_items ~ as.factor(trt) + lr_1 + unemployed + as.factor(trt) * sens, data=data3_age4, se_type="HC3")

idt <- c("Identity Threat\nvs. Control", lm1i$coefficients[7]) 
idt[2] <- lapply(idt[2], as.numeric); idt <- as.data.frame(idt)
colnames(idt) <- c("Treatment_Type", "Effect")
idt$se <- (lm1i$std.error)[7]
idt <- idt %>% mutate(cih90 = Effect + 1.645*se)
idt <- idt %>% mutate(cil90 = Effect - 1.645*se)
idt <- idt %>% mutate(cih95 = Effect + 1.96*se)
idt <- idt %>% mutate(cil95 = Effect - 1.96*se)

cit <- c("Whataboutism\nvs. Identity Threat", lm2i$coefficients[8]) 
cit[2] <- lapply(cit[2], as.numeric); cit <- as.data.frame(cit)
colnames(cit) <- c("Treatment_Type", "Effect")
cit$se <- (lm2i$std.error)[8]
cit <- cit %>% mutate(cih90 = Effect + 1.645*se)
cit <- cit %>% mutate(cil90 = Effect - 1.645*se)
cit <- cit %>% mutate(cih95 = Effect + 1.96*se)
cit <- cit %>% mutate(cil95 = Effect - 1.96*se)

com <- c("Whataboutism\nvs. Control", lm3i$coefficients[6]) 
com[2] <- lapply(com[2], as.numeric); com <- as.data.frame(com)
colnames(com) <- c("Treatment Type", "Effect")
com$se <- (lm3i$std.error)[6]
com <- com %>% mutate(cih90 = Effect + 1.645*se)
com <- com %>% mutate(cil90 = Effect - 1.645*se)
com <- com %>% mutate(cih95 = Effect + 1.96*se)
com <- com %>% mutate(cil95 = Effect - 1.96*se)

#bind results by treatment type
results_age4 <- rbind(unlist(idt), unlist(cit), unlist(com))
colnames(results_age4) <- c("Treatment_Type", "Effect", "se", "cih90", "cil90", "cih95", "cil95")
results_age4 <- data.frame(results_age4)
results_age4$Effect <- as.numeric(results_age4$Effect)
results_age4$cil90 <- as.numeric(results_age4$cil90)
results_age4$cih90 <- as.numeric(results_age4$cih90)
results_age4$cil95 <- as.numeric(results_age4$cil95)
results_age4$cih95 <- as.numeric(results_age4$cih95)
results_age4$Age<- "51-60"

#age group 5

lm1i <- lm_robust(list_items ~ as.factor(trt) + education + income + newjob + as.factor(trt) * sens, data=data1_age5, se_type="HC3")
lm2i <- lm_robust(list_items ~ as.factor(trt) + lr_1 + natid_covar1_1 + income + newjob + as.factor(trt) * sens, data=data2_age5, se_type="HC3")
lm3i <- lm_robust(list_items ~ as.factor(trt) + lr_1 + unemployed + as.factor(trt) * sens, data=data3_age5, se_type="HC3")

idt <- c("Identity Threat\nvs. Control", lm1i$coefficients[7]) 
idt[2] <- lapply(idt[2], as.numeric); idt <- as.data.frame(idt)
colnames(idt) <- c("Treatment_Type", "Effect")
idt$se <- (lm1i$std.error)[7]
idt <- idt %>% mutate(cih90 = Effect + 1.645*se)
idt <- idt %>% mutate(cil90 = Effect - 1.645*se)
idt <- idt %>% mutate(cih95 = Effect + 1.96*se)
idt <- idt %>% mutate(cil95 = Effect - 1.96*se)

cit <- c("Whataboutism\nvs. Identity Threat", lm2i$coefficients[8]) 
cit[2] <- lapply(cit[2], as.numeric); cit <- as.data.frame(cit)
colnames(cit) <- c("Treatment_Type", "Effect")
cit$se <- (lm2i$std.error)[8]
cit <- cit %>% mutate(cih90 = Effect + 1.645*se)
cit <- cit %>% mutate(cil90 = Effect - 1.645*se)
cit <- cit %>% mutate(cih95 = Effect + 1.96*se)
cit <- cit %>% mutate(cil95 = Effect - 1.96*se)

com <- c("Whataboutism\nvs. Control", lm3i$coefficients[6]) 
com[2] <- lapply(com[2], as.numeric); com <- as.data.frame(com)
colnames(com) <- c("Treatment Type", "Effect")
com$se <- (lm3i$std.error)[6]
com <- com %>% mutate(cih90 = Effect + 1.645*se)
com <- com %>% mutate(cil90 = Effect - 1.645*se)
com <- com %>% mutate(cih95 = Effect + 1.96*se)
com <- com %>% mutate(cil95 = Effect - 1.96*se)

#bind results by treatment type
results_age5 <- rbind(unlist(idt), unlist(cit), unlist(com))
colnames(results_age5) <- c("Treatment_Type", "Effect", "se", "cih90", "cil90", "cih95", "cil95")
results_age5 <- data.frame(results_age5)
results_age5$Effect <- as.numeric(results_age5$Effect)
results_age5$cil90 <- as.numeric(results_age5$cil90)
results_age5$cih90 <- as.numeric(results_age5$cih90)
results_age5$cil95 <- as.numeric(results_age5$cil95)
results_age5$cih95 <- as.numeric(results_age5$cih95)
results_age5$Age<- "61+"

#combine age groups and plot
age_hte <- rbind(results_age1, results_age2, results_age3, results_age4, results_age5)
age_hte$Age <- factor(age_hte$Age, levels = c("18-30", "31-40", "41-50", "51-60", "61+"))

ip2 <- ggplot(age_hte, aes(x = Age, y = Effect, colour=factor(Treatment_Type, level=factor_levels))) +
  geom_point(stat="identity", position=position_dodge(width=0.5), size=3) + 
  labs(tag = "B") +
  scale_y_continuous(limits = c(-0.8, 0.8)) +
  geom_linerange(aes(ymin = cil95, ymax = cih95), 
                 position = position_dodge(width=0.5,),
                 size = 1) +
  geom_linerange(aes(ymin = cil90, ymax = cih90), 
                 position = position_dodge(width=0.5),
                 size = 2) +
  ylab("Treatment Effect") +
  xlab("Age") +
  geom_hline(yintercept=0, linetype="dashed", color="blue") +
  labs(color='Update Type') +
  theme_bw()+
  theme(
    axis.text=element_text(size=16),
    plot.title = element_text(size=25),
    axis.title=element_text(size=18),
    legend.text=element_text(size=12),
    legend.title=element_text(size=16),
    legend.key.height = unit(1.5, 'cm')
  )

####TABLE D2(C) - Implicit Prejudice by Gender####

#women

#models by treatment type
lm1i <- lm_robust(list_items ~ as.factor(trt) + education + income + newjob + as.factor(trt) * sens, data = data1_f, se_type="HC3")
lm2i <- lm_robust(list_items ~ as.factor(trt) + lr_1 + natid_covar1_1 + income + newjob + as.factor(trt) * sens, data = data2_f, se_type="HC3")
lm3i <- lm_robust(list_items ~ as.factor(trt) + lr_1 + unemployed + as.factor(trt) * sens, data = data3_f, se_type="HC3")

# calculate confidence intervals
com1 <- c(
  "Identity Threat\nvs. Control", 
  lm1i$coefficients[7], 
  lm1i$coefficients[7] + z_95 * lm1i$std.error[7], 
  lm1i$coefficients[7] - z_95 * lm1i$std.error[7], 
  lm1i$coefficients[7] + z_90 * lm1i$std.error[7], 
  lm1i$coefficients[7] - z_90 * lm1i$std.error[7]
)
com2 <- c(
  "Whataboutism\nvs. Identity Threat", 
  lm2i$coefficients[8], 
  lm2i$coefficients[8] + z_95 * lm2i$std.error[8], 
  lm2i$coefficients[8] - z_95 * lm2i$std.error[8], 
  lm2i$coefficients[8] + z_90 * lm2i$std.error[8], 
  lm2i$coefficients[8] - z_90 * lm2i$std.error[8]
)
com3 <- c(
  "Whataboutism\nvs. Control", 
  lm3i$coefficients[6], 
  lm3i$coefficients[6] + z_95 * lm3i$std.error[6], 
  lm3i$coefficients[6] - z_95 * lm3i$std.error[6], 
  lm3i$coefficients[6] + z_90 * lm3i$std.error[6], 
  lm3i$coefficients[6] - z_90 * lm3i$std.error[6]
)

#bind results by treatment type
com_f <- data.frame(rbind(com1, com2, com3))
com_f$Gender <- "Female"

#male

#models by treatment type
lm1i <- lm_robust(list_items ~ as.factor(trt) + education + income + newjob + as.factor(trt) * sens, data = data1_m, se_type="HC3")
lm2i <- lm_robust(list_items ~ as.factor(trt) + lr_1 + natid_covar1_1 + income + newjob + as.factor(trt) * sens, data = data2_m, se_type="HC3")
lm3i <- lm_robust(list_items ~ as.factor(trt) + lr_1 + unemployed + as.factor(trt) * sens, data = data3_m, se_type="HC3")

# calculate confidence intervals
com1 <- c(
  "Identity Threat\nvs. Control", 
  lm1i$coefficients[7], 
  lm1i$coefficients[7] + z_95 * lm1i$std.error[7], 
  lm1i$coefficients[7] - z_95 * lm1i$std.error[7], 
  lm1i$coefficients[7] + z_90 * lm1i$std.error[7], 
  lm1i$coefficients[7] - z_90 * lm1i$std.error[7]
)
com2 <- c(
  "Whataboutism\nvs. Identity Threat", 
  lm2i$coefficients[8], 
  lm2i$coefficients[8] + z_95 * lm2i$std.error[8], 
  lm2i$coefficients[8] - z_95 * lm2i$std.error[8], 
  lm2i$coefficients[8] + z_90 * lm2i$std.error[8], 
  lm2i$coefficients[8] - z_90 * lm2i$std.error[8]
)
com3 <- c(
  "Whataboutism\nvs. Control", 
  lm3i$coefficients[6], 
  lm3i$coefficients[6] + z_95 * lm3i$std.error[6], 
  lm3i$coefficients[6] - z_95 * lm3i$std.error[6], 
  lm3i$coefficients[6] + z_90 * lm3i$std.error[6], 
  lm3i$coefficients[6] - z_90 * lm3i$std.error[6]
)

#bind results by treatment type
com_m <- data.frame(rbind(com1, com2, com3))
com_m$Gender <- "Male"

# Combine female/male and plot
gender_hte <- rbind(com_f, com_m)
colnames(gender_hte) <- c("treatment", "effect", "cih95", "cil95", "cih90", "cil90", "Gender")

# Ensure numeric data types for relevant columns
gender_hte$effect <- as.numeric(gender_hte$effect)
gender_hte$cih95 <- as.numeric(gender_hte$cih95)
gender_hte$cil95 <- as.numeric(gender_hte$cil95)
gender_hte$cih90 <- as.numeric(gender_hte$cih90)
gender_hte$cil90 <- as.numeric(gender_hte$cil90)

ip3 <- ggplot(gender_hte, aes(x = Gender, y = effect, colour=factor(treatment, level=factor_levels))) +
  geom_point(stat="identity", position=position_dodge(width=0.3), size=3) +
  labs(tag = "C") +
  scale_y_continuous(limits = c(-0.8, 0.8)) +
  geom_linerange(aes(ymin = cil95, ymax = cih95), 
                 position = position_dodge(width=0.3,),
                 size = 1) +
  geom_linerange(aes(ymin = cil90, ymax = cih90), 
                 position = position_dodge(width=0.3),
                 size = 2) +
  ylab("Treatment Effect") +
  xlab("Gender") +
  geom_hline(yintercept=0, linetype="dashed", color="blue") +
  labs(color='Update Type') +
  theme_bw()+
  theme(
    axis.text=element_text(size=16),
    plot.title = element_text(size=25),
    axis.title=element_text(size=18),
    legend.text=element_text(size=12),
    legend.title=element_text(size=16),
    legend.key.height = unit(1.5, 'cm')
  )

####TABLE D2(D) - Implicit Prejudice by Pre-Treatment National Identification####

#quintile 1

#models by treatment type
lm1i <- lm_robust(list_items ~ as.factor(trt) + education + income + newjob + as.factor(trt) * sens, data = data1_nid1, se_type="HC3")
lm2i <- lm_robust(list_items ~ as.factor(trt) + lr_1 + natid_covar1_1 + income + newjob + as.factor(trt) * sens, data = data2_nid1, se_type="HC3")
lm3i <- lm_robust(list_items ~ as.factor(trt) + lr_1 + unemployed + as.factor(trt) * sens, data = data3_nid1, se_type="HC3")

idt <- c("Identity Threat\nvs. Control", lm1i$coefficients[7]) 
idt[2] <- lapply(idt[2], as.numeric); idt <- as.data.frame(idt)
colnames(idt) <- c("Treatment_Type", "Effect")
idt$se <- (lm1i$std.error)[7]
idt <- idt %>% mutate(cih95 = Effect + 1.96*se)
idt <- idt %>% mutate(cil95 = Effect - 1.96*se)
idt <- idt %>% mutate(cih90 = Effect + 1.645*se)
idt <- idt %>% mutate(cil90 = Effect - 1.645*se)

cit <- c("Whataboutism\nvs. Identity Threat", lm2i$coefficients[8]) 
cit[2] <- lapply(cit[2], as.numeric); cit <- as.data.frame(cit)
colnames(cit) <- c("Treatment_Type", "Effect")
cit$se <- (lm2i$std.error)[8]
cit <- cit %>% mutate(cih95 = Effect + 1.96*se)
cit <- cit %>% mutate(cil95 = Effect - 1.96*se)
cit <- cit %>% mutate(cih90 = Effect + 1.645*se)
cit <- cit %>% mutate(cil90 = Effect - 1.645*se)

com <- c("Whataboutism\nvs. Control", lm3i$coefficients[6]) 
com[2] <- lapply(com[2], as.numeric); com <- as.data.frame(com)
colnames(com) <- c("Treatment Type", "Effect")
com$se <- (lm3i$std.error)[6]
com <- com %>% mutate(cih95 = Effect + 1.96*se)
com <- com %>% mutate(cil95 = Effect - 1.96*se)
com <- com %>% mutate(cih90 = Effect + 1.645*se)
com <- com %>% mutate(cil90 = Effect - 1.645*se)

#bind results by treatment type
results_nid1 <- rbind(unlist(idt), unlist(cit), unlist(com))
colnames(results_nid1) <- c("Treatment_Type", "Effect", "se", "cih95", "cil95", "cih90", "cil90")
results_nid1 <- data.frame(results_nid1)
results_nid1$Effect <- as.numeric(results_nid1$Effect)
results_nid1$cil95 <- as.numeric(results_nid1$cil95)
results_nid1$cil90 <- as.numeric(results_nid1$cil90)
results_nid1$cih95 <- as.numeric(results_nid1$cih95)
results_nid1$cih90 <- as.numeric(results_nid1$cih90)
results_nid1$Quintile<- 1

#quintile 2

#models by treatment type
lm1i <- lm_robust(list_items ~ as.factor(trt) + education + income + newjob + as.factor(trt) * sens, data = data1_nid2, se_type="HC3")
lm2i <- lm_robust(list_items ~ as.factor(trt) + lr_1 + natid_covar1_1 + income + newjob + as.factor(trt) * sens, data = data2_nid2, se_type="HC3")
lm3i <- lm_robust(list_items ~ as.factor(trt) + lr_1 + unemployed + as.factor(trt) * sens, data = data3_nid2, se_type="HC3")

idt <- c("Identity Threat\nvs. Control", lm1i$coefficients[7]) 
idt[2] <- lapply(idt[2], as.numeric); idt <- as.data.frame(idt)
colnames(idt) <- c("Treatment_Type", "Effect")
idt$se <- (lm1i$std.error)[7]
idt <- idt %>% mutate(cih95 = Effect + 1.96*se)
idt <- idt %>% mutate(cil95 = Effect - 1.96*se)
idt <- idt %>% mutate(cih90 = Effect + 1.645*se)
idt <- idt %>% mutate(cil90 = Effect - 1.645*se)

cit <- c("Whataboutism\nvs. Identity Threat", lm2i$coefficients[8]) 
cit[2] <- lapply(cit[2], as.numeric); cit <- as.data.frame(cit)
colnames(cit) <- c("Treatment_Type", "Effect")
cit$se <- (lm2i$std.error)[8]
cit <- cit %>% mutate(cih95 = Effect + 1.96*se)
cit <- cit %>% mutate(cil95 = Effect - 1.96*se)
cit <- cit %>% mutate(cih90 = Effect + 1.645*se)
cit <- cit %>% mutate(cil90 = Effect - 1.645*se)

com <- c("Whataboutism\nvs. Control", lm3i$coefficients[6]) 
com[2] <- lapply(com[2], as.numeric); com <- as.data.frame(com)
colnames(com) <- c("Treatment_Type", "Effect")
com$se <- (lm3i$std.error)[6]
com <- com %>% mutate(cih95 = Effect + 1.96*se)
com <- com %>% mutate(cil95 = Effect - 1.96*se)
com <- com %>% mutate(cih90 = Effect + 1.645*se)
com <- com %>% mutate(cil90 = Effect - 1.645*se)

#bind results by treatment type
results_nid2 <- rbind(unlist(idt), unlist(cit), unlist(com))
colnames(results_nid2) <- c("Treatment_Type", "Effect", "se", "cih95", "cil95", "cih90", "cil90")
results_nid2 <- data.frame(results_nid2)
results_nid2$Effect <- as.numeric(results_nid2$Effect)
results_nid2$cil95 <- as.numeric(results_nid2$cil95)
results_nid2$cil90 <- as.numeric(results_nid2$cil90)
results_nid2$cih95 <- as.numeric(results_nid2$cih95)
results_nid2$cih90 <- as.numeric(results_nid2$cih90)
results_nid2$Quintile<- 2

#quintile 3

#models by treatment type
lm1i <- lm_robust(list_items ~ as.factor(trt) + education + income + newjob + as.factor(trt) * sens, data = data1_nid3, se_type="HC3")
lm2i <- lm_robust(list_items ~ as.factor(trt) + lr_1 + natid_covar1_1 + income + newjob + as.factor(trt) * sens, data = data2_nid3, se_type="HC3")
lm3i <- lm_robust(list_items ~ as.factor(trt) + lr_1 + unemployed + as.factor(trt) * sens, data = data3_nid3, se_type="HC3")

idt <- c("Identity Threat\nvs. Control", lm1i$coefficients[7]) 
idt[2] <- lapply(idt[2], as.numeric); idt <- as.data.frame(idt)
colnames(idt) <- c("Treatment_Type", "Effect")
idt$se <- (lm1i$std.error)[7]
idt <- idt %>% mutate(cih95 = Effect + 1.96*se)
idt <- idt %>% mutate(cil95 = Effect - 1.96*se)
idt <- idt %>% mutate(cih90 = Effect + 1.645*se)
idt <- idt %>% mutate(cil90 = Effect - 1.645*se)

cit <- c("Whataboutism\nvs. Identity Threat", lm2i$coefficients[8]) 
cit[2] <- lapply(cit[2], as.numeric); cit <- as.data.frame(cit)
colnames(cit) <- c("Treatment_Type", "Effect")
cit$se <- (lm2i$std.error)[8]
cit <- cit %>% mutate(cih95 = Effect + 1.96*se)
cit <- cit %>% mutate(cil95 = Effect - 1.96*se)
cit <- cit %>% mutate(cih90 = Effect + 1.645*se)
cit <- cit %>% mutate(cil90 = Effect - 1.645*se)

com <- c("Whataboutism\nvs. Control", lm3i$coefficients[6]) 
com[2] <- lapply(com[2], as.numeric); com <- as.data.frame(com)
colnames(com) <- c("Treatment_Type", "Effect")
com$se <- (lm3i$std.error)[6]
com <- com %>% mutate(cih95 = Effect + 1.96*se)
com <- com %>% mutate(cil95 = Effect - 1.96*se)
com <- com %>% mutate(cih90 = Effect + 1.645*se)
com <- com %>% mutate(cil90 = Effect - 1.645*se)

#bind results by treatment type
results_nid3 <- rbind(unlist(idt), unlist(cit), unlist(com))
colnames(results_nid3) <- c("Treatment_Type", "Effect", "se", "cih95", "cil95", "cih90", "cil90")
results_nid3 <- data.frame(results_nid3)
results_nid3$Effect <- as.numeric(results_nid3$Effect)
results_nid3$cil95 <- as.numeric(results_nid3$cil95)
results_nid3$cil90 <- as.numeric(results_nid3$cil90)
results_nid3$cih95 <- as.numeric(results_nid3$cih95)
results_nid3$cih90 <- as.numeric(results_nid3$cih90)
results_nid3$Quintile<- 3

#quintile 4

#models by treatment type
lm1i <- lm_robust(list_items ~ as.factor(trt) + education + income + newjob + as.factor(trt) * sens, data = data1_nid4, se_type="HC3")
lm2i <- lm_robust(list_items ~ as.factor(trt) + lr_1 + natid_covar1_1 + income + newjob + as.factor(trt) * sens, data = data2_nid4, se_type="HC3")
lm3i <- lm_robust(list_items ~ as.factor(trt) + lr_1 + unemployed + as.factor(trt) * sens, data = data3_nid4, se_type="HC3")

idt <- c("Identity Threat\nvs. Control", lm1i$coefficients[7]) 
idt[2] <- lapply(idt[2], as.numeric); idt <- as.data.frame(idt)
colnames(idt) <- c("Treatment_Type", "Effect")
idt$se <- (lm1i$std.error)[7]
idt <- idt %>% mutate(cih95 = Effect + 1.96*se)
idt <- idt %>% mutate(cil95 = Effect - 1.96*se)
idt <- idt %>% mutate(cih90 = Effect + 1.645*se)
idt <- idt %>% mutate(cil90 = Effect - 1.645*se)

cit <- c("Whataboutism\nvs. Identity Threat", lm2i$coefficients[8]) 
cit[2] <- lapply(cit[2], as.numeric); cit <- as.data.frame(cit)
colnames(cit) <- c("Treatment_Type", "Effect")
cit$se <- (lm2i$std.error)[8]
cit <- cit %>% mutate(cih95 = Effect + 1.96*se)
cit <- cit %>% mutate(cil95 = Effect - 1.96*se)
cit <- cit %>% mutate(cih90 = Effect + 1.645*se)
cit <- cit %>% mutate(cil90 = Effect - 1.645*se)

com <- c("Whataboutism\nvs. Control", lm3i$coefficients[6]) 
com[2] <- lapply(com[2], as.numeric); com <- as.data.frame(com)
colnames(com) <- c("Treatment_Type", "Effect")
com$se <- (lm3i$std.error)[6]
com <- com %>% mutate(cih95 = Effect + 1.96*se)
com <- com %>% mutate(cil95 = Effect - 1.96*se)
com <- com %>% mutate(cih90 = Effect + 1.645*se)
com <- com %>% mutate(cil90 = Effect - 1.645*se)

#bind results by treatment type
results_nid4 <- rbind(unlist(idt), unlist(cit), unlist(com))
colnames(results_nid4) <- c("Treatment_Type", "Effect", "se", "cih95", "cil95", "cih90", "cil90")
results_nid4 <- data.frame(results_nid4)
results_nid4$Effect <- as.numeric(results_nid4$Effect)
results_nid4$cil95 <- as.numeric(results_nid4$cil95)
results_nid4$cil90 <- as.numeric(results_nid4$cil90)
results_nid4$cih95 <- as.numeric(results_nid4$cih95)
results_nid4$cih90 <- as.numeric(results_nid4$cih90)
results_nid4$Quintile<- 4

#quintile 5

#models by treatment type
lm1i <- lm_robust(list_items ~ as.factor(trt) + education + income + newjob + as.factor(trt) * sens, data = data1_nid5, se_type="HC3")
lm2i <- lm_robust(list_items ~ as.factor(trt) + lr_1 + natid_covar1_1 + income + newjob + as.factor(trt) * sens, data = data2_nid5, se_type="HC3")
lm3i <- lm_robust(list_items ~ as.factor(trt) + lr_1 + unemployed + as.factor(trt) * sens, data = data3_nid5, se_type="HC3")

idt <- c("Identity Threat\nvs. Control", lm1i$coefficients[7]) 
idt[2] <- lapply(idt[2], as.numeric); idt <- as.data.frame(idt)
colnames(idt) <- c("Treatment_Type", "Effect")
idt$se <- (lm1i$std.error)[7]
idt <- idt %>% mutate(cih95 = Effect + 1.96*se)
idt <- idt %>% mutate(cil95 = Effect - 1.96*se)
idt <- idt %>% mutate(cih90 = Effect + 1.645*se)
idt <- idt %>% mutate(cil90 = Effect - 1.645*se)

cit <- c("Whataboutism\nvs. Identity Threat", lm2i$coefficients[8]) 
cit[2] <- lapply(cit[2], as.numeric); cit <- as.data.frame(cit)
colnames(cit) <- c("Treatment_Type", "Effect")
cit$se <- (lm2i$std.error)[8]
cit <- cit %>% mutate(cih95 = Effect + 1.96*se)
cit <- cit %>% mutate(cil95 = Effect - 1.96*se)
cit <- cit %>% mutate(cih90 = Effect + 1.645*se)
cit <- cit %>% mutate(cil90 = Effect - 1.645*se)

com <- c("Whataboutism\nvs. Control", lm3i$coefficients[6]) 
com[2] <- lapply(com[2], as.numeric); com <- as.data.frame(com)
colnames(com) <- c("Treatment_Type", "Effect")
com$se <- (lm3i$std.error)[6]
com <- com %>% mutate(cih95 = Effect + 1.96*se)
com <- com %>% mutate(cil95 = Effect - 1.96*se)
com <- com %>% mutate(cih90 = Effect + 1.645*se)
com <- com %>% mutate(cil90 = Effect - 1.645*se)

#bind results by treatment type
results_nid5 <- rbind(unlist(idt), unlist(cit), unlist(com))
colnames(results_nid5) <- c("Treatment_Type", "Effect", "se", "cih95", "cil95", "cih90", "cil90")
results_nid5 <- data.frame(results_nid5)
results_nid5$Effect <- as.numeric(results_nid5$Effect)
results_nid5$cil95 <- as.numeric(results_nid5$cil95)
results_nid5$cil90 <- as.numeric(results_nid5$cil90)
results_nid5$cih95 <- as.numeric(results_nid5$cih95)
results_nid5$cih90 <- as.numeric(results_nid5$cih90)
results_nid5$Quintile<- 5

#combine quintiles and plot
nid_hte <- rbind(results_nid1, results_nid2, results_nid3, results_nid4, results_nid5)
nid_hte$Effect <- as.numeric(as.character(nid_hte$Effect))

ip4 <- ggplot(nid_hte, aes(x = Quintile, y = Effect, colour=factor(Treatment_Type, level=factor_levels))) +
  geom_point(stat="identity", position=position_dodge(width=0.5), size=3) + 
  labs(tag = "D") +
  scale_y_continuous(limits = c(-0.8, 0.8)) +
  geom_linerange(aes(ymin = cil95, ymax = cih95), 
                 position = position_dodge(width=0.5,),
                 size = 1) +
  geom_linerange(aes(ymin = cil90, ymax = cih90), 
                 position = position_dodge(width=0.5),
                 size = 2)+
  ylab("Treatment Effect") +
  xlab("National Identification") +
  geom_hline(yintercept=0, linetype="dashed", color="blue") +
  labs(color='Update Type') +
  theme_bw()+
  theme(
    axis.text=element_text(size=16),
    plot.title = element_text(size=25),
    axis.title=element_text(size=18),
    legend.text=element_text(size=12),
    legend.title=element_text(size=16),
    legend.key.height = unit(1.5, 'cm')
  )

####Combined Panel for Figure D2####

(ip1 | ip2) /
  (ip3 | ip4)+
  plot_layout(guides = "collect") &
  theme(legend.position = "right")

ggsave("figD2.png", width = 12, height = 10, dpi = 300)

####TABLE D3 - HTEs for BLM Donation####
####TABLE D3(A) - BLM Donation by Left/Right Self-Placement####

#Left-wing
data1_left <- subset(data1, data1$lr_1 < 5)
data2_left <- subset(data2, data2$lr_1 < 5)
data3_left <- subset(data3, data3$lr_1 < 5)

#models by treatment type
lm1a <- lm_robust(donation_amount ~ as.factor(trt) + education + income + newjob, data=data1_left, se_type="HC3")
lm2a <- lm_robust(donation_amount ~ as.factor(trt) + lr_1 + natid_covar1_1 + income + newjob, data = data2_left, se_type="HC3")
lm3a <- lm_robust(donation_amount ~ as.factor(trt) + lr_1 + unemployed, data = data3_left, se_type="HC3")

idt <- c("Identity Threat\nvs. Control", lm1a$coefficients[2])
idt[2] <- lapply(idt[2], as.numeric); idt <- as.data.frame(idt)
colnames(idt) <- c("Treatment_Type", "Effect")
idt$se <- lm1a$std.error[2]
idt <- idt %>% 
  mutate(
    cih_95 = Effect + z_95 * se,
    cil_95 = Effect - z_95 * se,
    cih_90 = Effect + z_90 * se,
    cil_90 = Effect - z_90 * se
  )

cit <- c("Whataboutism\nvs. Identity Threat", lm2a$coefficients[2])
cit[2] <- lapply(cit[2], as.numeric); cit <- as.data.frame(cit)
colnames(cit) <- c("Treatment_Type", "Effect")
cit$se <- lm2a$std.error[2]
cit <- cit %>% 
  mutate(
    cih_95 = Effect + z_95 * se,
    cil_95 = Effect - z_95 * se,
    cih_90 = Effect + z_90 * se,
    cil_90 = Effect - z_90 * se
  )

com <- c("Whataboutism\nvs. Control", lm3a$coefficients[2])
com[2] <- lapply(com[2], as.numeric); com <- as.data.frame(com)
colnames(com) <- c("Treatment_Type", "Effect")
com$se <- lm3a$std.error[2]
com <- com %>% 
  mutate(
    cih_95 = Effect + z_95 * se,
    cil_95 = Effect - z_95 * se,
    cih_90 = Effect + z_90 * se,
    cil_90 = Effect - z_90 * se
  )

#bind results by treatment type
results_left <- rbind(unlist(idt), unlist(cit), unlist(com))
colnames(results_left) <- c("Treatment_Type", "Effect", "se", "cih_95", "cil_95", "cih_90", "cil_90")
results_left <- data.frame(results_left)
results_left$Effect <- as.numeric(results_left$Effect)
results_left$cil_95 <- as.numeric(results_left$cil_95)
results_left$cih_95 <- as.numeric(results_left$cih_95)
results_left$cil_90 <- as.numeric(results_left$cil_90)
results_left$cih_90 <- as.numeric(results_left$cih_90)
results_left$HTE <- "Left"

#Right-wing
data1_right <- subset(data1, data1$lr_1 >= 6)
data2_right <- subset(data2, data2$lr_1 >= 6)
data3_right <- subset(data3, data3$lr_1 >= 6)

#models by treatment type
lm1a <- lm_robust(donation_amount ~ as.factor(trt) + education + income + newjob, data=data1_right, se_type="HC3")
lm2a <- lm_robust(donation_amount ~ as.factor(trt) + lr_1 + natid_covar1_1 + income + newjob, data = data2_right, se_type="HC3")
lm3a <- lm_robust(donation_amount ~ as.factor(trt) + lr_1 + unemployed, data = data3_right, se_type="HC3")

idt <- c("Identity Threat\nvs. Control", lm1a$coefficients[2])
idt[2] <- lapply(idt[2], as.numeric); idt <- as.data.frame(idt)
colnames(idt) <- c("Treatment_Type", "Effect")
idt$se <- lm1a$std.error[2]
idt <- idt %>% 
  mutate(
    cih_95 = Effect + z_95 * se,
    cil_95 = Effect - z_95 * se,
    cih_90 = Effect + z_90 * se,
    cil_90 = Effect - z_90 * se
  )

cit <- c("Whataboutism\nvs. Identity Threat", lm2a$coefficients[2])
cit[2] <- lapply(cit[2], as.numeric); cit <- as.data.frame(cit)
colnames(cit) <- c("Treatment_Type", "Effect")
cit$se <- lm2a$std.error[2]
cit <- cit %>% 
  mutate(
    cih_95 = Effect + z_95 * se,
    cil_95 = Effect - z_95 * se,
    cih_90 = Effect + z_90 * se,
    cil_90 = Effect - z_90 * se
  )

com <- c("Whataboutism\nvs. Control", lm3a$coefficients[2])
com[2] <- lapply(com[2], as.numeric); com <- as.data.frame(com)
colnames(com) <- c("Treatment_Type", "Effect")
com$se <- lm3a$std.error[2]
com <- com %>% 
  mutate(
    cih_95 = Effect + z_95 * se,
    cil_95 = Effect - z_95 * se,
    cih_90 = Effect + z_90 * se,
    cil_90 = Effect - z_90 * se
  )

#bind results by treatment type
results_right <- rbind(unlist(idt), unlist(cit), unlist(com))
colnames(results_right) <- c("Treatment_Type", "Effect", "se", "cih_95", "cil_95", "cih_90", "cil_90")
results_right <- data.frame(results_right)
results_right$Effect <- as.numeric(results_right$Effect)
results_right$cil_95 <- as.numeric(results_right$cil_95)
results_right$cih_95 <- as.numeric(results_right$cih_95)
results_right$cil_90 <- as.numeric(results_right$cil_90)
results_right$cih_90 <- as.numeric(results_right$cih_90)
results_right$HTE <- "Right"

#bind left and right results and plot
lr_hte <- rbind(results_left, results_right)
lr_hte$HTE <- factor(lr_hte$HTE, levels = c("Left", "Right"))

do1 <- ggplot(lr_hte, aes(x = HTE, y = Effect, colour=factor(Treatment_Type, level=factor_levels))) +
  geom_point(stat="identity", position=position_dodge(width=0.3), size=3) + 
  labs(tag = "A") +
  scale_y_continuous(limits = c(-7.2, 6.6)) +
  geom_linerange(aes(ymin = cil_95, ymax = cih_95), 
                 position = position_dodge(width=0.3,),
                 size = 1) +
  geom_linerange(aes(ymin = cil_90, ymax = cih_90), 
                 position = position_dodge(width=0.3),
                 size = 2)+
  ylab("Treatment Effect") +
  xlab("Left-Right Self-Placement") +
  geom_hline(yintercept=0, linetype="dashed", color="blue") +
  labs(color='Update Type') +
  theme_bw()+
  theme(
    axis.text=element_text(size=16),
    plot.title = element_text(size=25),
    axis.title=element_text(size=18),
    legend.text=element_text(size=12),
    legend.title=element_text(size=16),
    legend.key.height = unit(1.5, 'cm')
  )

####TABLE D3(B) - BLM Donation by Age####

#age group 1

data1_age1 <- subset(data1, data1$age<31)
data2_age1 <- subset(data2, data2$age<31)
data3_age1 <- subset(data3, data3$age<31)

#models by treatment type
lm1a <- lm_robust(donation_amount ~ as.factor(trt) + education + income + newjob, data=data1_age1, se_type="HC3")
lm2a <- lm_robust(donation_amount ~ as.factor(trt) + lr_1  + natid_covar1_1 + income + newjob, data=data2_age1, se_type="HC3")
lm3a <- lm_robust(donation_amount ~ as.factor(trt) + lr_1 + unemployed, data=data3_age1, se_type="HC3")

idt <- c("Identity Threat\nvs. Control", lm1a$coefficients[2]) 
idt[2] <- lapply(idt[2], as.numeric); idt <- as.data.frame(idt)
colnames(idt) <- c("Treatment_Type", "Effect")
idt$se <- (lm1a$std.error)[2]
idt <- idt %>% mutate(cih90 = Effect + 1.645*se)
idt <- idt %>% mutate(cil90 = Effect - 1.645*se)
idt <- idt %>% mutate(cih95 = Effect + 1.96*se)
idt <- idt %>% mutate(cil95 = Effect - 1.96*se)

cit <- c("Whataboutism\nvs. Identity Threat", lm2a$coefficients[2]) 
cit[2] <- lapply(cit[2], as.numeric); cit <- as.data.frame(cit)
colnames(cit) <- c("Treatment_Type", "Effect")
cit$se <- (lm2a$std.error)[2]
cit <- cit %>% mutate(cih90 = Effect + 1.645*se)
cit <- cit %>% mutate(cil90 = Effect - 1.645*se)
cit <- cit %>% mutate(cih95 = Effect + 1.96*se)
cit <- cit %>% mutate(cil95 = Effect - 1.96*se)

com <- c("Whataboutism\nvs. Control", lm3a$coefficients[2]) 
com[2] <- lapply(com[2], as.numeric); com <- as.data.frame(com)
colnames(com) <- c("Treatment_Type", "Effect")
com$se <- (lm3a$std.error)[2]
com <- com %>% mutate(cih90 = Effect + 1.645*se)
com <- com %>% mutate(cil90 = Effect - 1.645*se)
com <- com %>% mutate(cih95 = Effect + 1.96*se)
com <- com %>% mutate(cil95 = Effect - 1.96*se)

#bind results by treatment type
results_age1 <- rbind(unlist(idt), unlist(cit), unlist(com))
colnames(results_age1) <- c("Treatment_Type", "Effect", "se", "cih90", "cil90", "cih95", "cil95")
results_age1 <- data.frame(results_age1)
results_age1$Effect <- as.numeric(results_age1$Effect)
results_age1$cil90 <- as.numeric(results_age1$cil90)
results_age1$cih90 <- as.numeric(results_age1$cih90)
results_age1$cil95 <- as.numeric(results_age1$cil95)
results_age1$cih95 <- as.numeric(results_age1$cih95)

results_age1$Age<- "18-30"

#age group 2

data1_age2 <- subset(data1, data1$age>30 & data1$age<41)
data2_age2 <- subset(data2, data2$age>30 & data2$age<41)
data3_age2 <- subset(data3, data3$age>30 & data3$age<41)

#models by treatment type
lm1a <- lm_robust(donation_amount ~ as.factor(trt) + education + income + newjob, data=data1_age2, se_type="HC3")
lm2a <- lm_robust(donation_amount ~ as.factor(trt) + lr_1  + natid_covar1_1 + income + newjob, data=data2_age2, se_type="HC3")
lm3a <- lm_robust(donation_amount ~ as.factor(trt) + lr_1 + unemployed, data=data3_age2, se_type="HC3")

idt <- c("Identity Threat\nvs. Control", lm1a$coefficients[2]) 
idt[2] <- lapply(idt[2], as.numeric); idt <- as.data.frame(idt)
colnames(idt) <- c("Treatment_Type", "Effect")
idt$se <- (lm1a$std.error)[2]
idt <- idt %>% mutate(cih90 = Effect + 1.645*se)
idt <- idt %>% mutate(cil90 = Effect - 1.645*se)
idt <- idt %>% mutate(cih95 = Effect + 1.96*se)
idt <- idt %>% mutate(cil95 = Effect - 1.96*se)

cit <- c("Whataboutism\nvs. Identity Threat", lm2a$coefficients[2]) 
cit[2] <- lapply(cit[2], as.numeric); cit <- as.data.frame(cit)
colnames(cit) <- c("Treatment_Type", "Effect")
cit$se <- (lm2a$std.error)[2]
cit <- cit %>% mutate(cih90 = Effect + 1.645*se)
cit <- cit %>% mutate(cil90 = Effect - 1.645*se)
cit <- cit %>% mutate(cih95 = Effect + 1.96*se)
cit <- cit %>% mutate(cil95 = Effect - 1.96*se)

com <- c("Whataboutism\nvs. Control", lm3a$coefficients[2]) 
com[2] <- lapply(com[2], as.numeric); com <- as.data.frame(com)
colnames(com) <- c("Treatment_Type", "Effect")
com$se <- (lm3a$std.error)[2]
com <- com %>% mutate(cih90 = Effect + 1.645*se)
com <- com %>% mutate(cil90 = Effect - 1.645*se)
com <- com %>% mutate(cih95 = Effect + 1.96*se)
com <- com %>% mutate(cil95 = Effect - 1.96*se)

#bind results by treatment type
results_age2 <- rbind(unlist(idt), unlist(cit), unlist(com))
colnames(results_age2) <- c("Treatment_Type", "Effect", "se", "cih90", "cil90", "cih95", "cil95")
results_age2 <- data.frame(results_age2)
results_age2$Effect <- as.numeric(results_age2$Effect)
results_age2$cil90 <- as.numeric(results_age2$cil90)
results_age2$cih90 <- as.numeric(results_age2$cih90)
results_age2$cil95 <- as.numeric(results_age2$cil95)
results_age2$cih95 <- as.numeric(results_age2$cih95)
results_age2$Age<- "31-40"

#age group 3

data1_age3 <- subset(data1, data1$age>40 & data1$age<51)
data2_age3 <- subset(data2, data2$age>40 & data2$age<51)
data3_age3 <- subset(data3, data3$age>40 & data3$age<51)

#models by treatment type
lm1a <- lm_robust(donation_amount ~ as.factor(trt) + education + income + newjob, data=data1_age3, se_type="HC3")
lm2a <- lm_robust(donation_amount ~ as.factor(trt) + lr_1  + natid_covar1_1 + income + newjob, data=data2_age3, se_type="HC3")
lm3a <- lm_robust(donation_amount ~ as.factor(trt) + lr_1 + unemployed, data=data3_age3, se_type="HC3")

idt <- c("Identity Threat\nvs. Control", lm1a$coefficients[2]) 
idt[2] <- lapply(idt[2], as.numeric); idt <- as.data.frame(idt)
colnames(idt) <- c("Treatment_Type", "Effect")
idt$se <- (lm1a$std.error)[2]
idt <- idt %>% mutate(cih90 = Effect + 1.645*se)
idt <- idt %>% mutate(cil90 = Effect - 1.645*se)
idt <- idt %>% mutate(cih95 = Effect + 1.96*se)
idt <- idt %>% mutate(cil95 = Effect - 1.96*se)

cit <- c("Whataboutism\nvs. Identity Threat", lm2a$coefficients[2]) 
cit[2] <- lapply(cit[2], as.numeric); cit <- as.data.frame(cit)
colnames(cit) <- c("Treatment_Type", "Effect")
cit$se <- (lm2a$std.error)[2]
cit <- cit %>% mutate(cih90 = Effect + 1.645*se)
cit <- cit %>% mutate(cil90 = Effect - 1.645*se)
cit <- cit %>% mutate(cih95 = Effect + 1.96*se)
cit <- cit %>% mutate(cil95 = Effect - 1.96*se)

com <- c("Whataboutism\nvs. Control", lm3a$coefficients[2]) 
com[2] <- lapply(com[2], as.numeric); com <- as.data.frame(com)
colnames(com) <- c("Treatment_Type", "Effect")
com$se <- (lm3a$std.error)[2]
com <- com %>% mutate(cih90 = Effect + 1.645*se)
com <- com %>% mutate(cil90 = Effect - 1.645*se)
com <- com %>% mutate(cih95 = Effect + 1.96*se)
com <- com %>% mutate(cil95 = Effect - 1.96*se)

#bind results by treatment type
results_age3 <- rbind(unlist(idt), unlist(cit), unlist(com))
colnames(results_age3) <- c("Treatment_Type", "Effect", "se", "cih90", "cil90", "cih95", "cil95")
results_age3 <- data.frame(results_age3)
results_age3$Effect <- as.numeric(results_age3$Effect)
results_age3$cil90 <- as.numeric(results_age3$cil90)
results_age3$cih90 <- as.numeric(results_age3$cih90)
results_age3$cil95 <- as.numeric(results_age3$cil95)
results_age3$cih95 <- as.numeric(results_age3$cih95)
results_age3$Age<- "41-50"

#age group 4

data1_age4 <- subset(data1, data1$age>50 & data1$age<61)
data2_age4 <- subset(data2, data2$age>50 & data2$age<61)
data3_age4 <- subset(data3, data3$age>50 & data3$age<61)

#models by treatment type
lm1a <- lm_robust(donation_amount ~ as.factor(trt) + education + income + newjob, data=data1_age4, se_type="HC3")
lm2a <- lm_robust(donation_amount ~ as.factor(trt) + lr_1  + natid_covar1_1 + income + newjob, data=data2_age4, se_type="HC3")
lm3a <- lm_robust(donation_amount ~ as.factor(trt) + lr_1 + unemployed, data=data3_age4, se_type="HC3")

idt <- c("Identity Threat\nvs. Control", lm1a$coefficients[2]) 
idt[2] <- lapply(idt[2], as.numeric); idt <- as.data.frame(idt)
colnames(idt) <- c("Treatment_Type", "Effect")
idt$se <- (lm1a$std.error)[2]
idt <- idt %>% mutate(cih90 = Effect + 1.645*se)
idt <- idt %>% mutate(cil90 = Effect - 1.645*se)
idt <- idt %>% mutate(cih95 = Effect + 1.96*se)
idt <- idt %>% mutate(cil95 = Effect - 1.96*se)

cit <- c("Whataboutism\nvs. Identity Threat", lm2a$coefficients[2]) 
cit[2] <- lapply(cit[2], as.numeric); cit <- as.data.frame(cit)
colnames(cit) <- c("Treatment_Type", "Effect")
cit$se <- (lm2a$std.error)[2]
cit <- cit %>% mutate(cih90 = Effect + 1.645*se)
cit <- cit %>% mutate(cil90 = Effect - 1.645*se)
cit <- cit %>% mutate(cih95 = Effect + 1.96*se)
cit <- cit %>% mutate(cil95 = Effect - 1.96*se)

com <- c("Whataboutism\nvs. Control", lm3a$coefficients[2]) 
com[2] <- lapply(com[2], as.numeric); com <- as.data.frame(com)
colnames(com) <- c("Treatment_Type", "Effect")
com$se <- (lm3a$std.error)[2]
com <- com %>% mutate(cih90 = Effect + 1.645*se)
com <- com %>% mutate(cil90 = Effect - 1.645*se)
com <- com %>% mutate(cih95 = Effect + 1.96*se)
com <- com %>% mutate(cil95 = Effect - 1.96*se)

#bind results by treatment type
results_age4 <- rbind(unlist(idt), unlist(cit), unlist(com))
colnames(results_age4) <- c("Treatment_Type", "Effect", "se", "cih90", "cil90", "cih95", "cil95")
results_age4 <- data.frame(results_age4)
results_age4$Effect <- as.numeric(results_age4$Effect)
results_age4$cil90 <- as.numeric(results_age4$cil90)
results_age4$cih90 <- as.numeric(results_age4$cih90)
results_age4$cil95 <- as.numeric(results_age4$cil95)
results_age4$cih95 <- as.numeric(results_age4$cih95)
results_age4$Age<- "51-60"

#age group 5

data1_age5 <- subset(data1, data1$age>60)
data2_age5 <- subset(data2, data2$age>60)
data3_age5 <- subset(data3, data3$age>60)

#models by treatment type
lm1a <- lm_robust(donation_amount ~ as.factor(trt) + education + income + newjob, data=data1_age5, se_type="HC3")
lm2a <- lm_robust(donation_amount ~ as.factor(trt) + lr_1  + natid_covar1_1 + income + newjob, data=data2_age5, se_type="HC3")
lm3a <- lm_robust(donation_amount ~ as.factor(trt) + lr_1 + unemployed, data=data3_age5, se_type="HC3")

idt <- c("Identity Threat\nvs. Control", lm1a$coefficients[2]) 
idt[2] <- lapply(idt[2], as.numeric); idt <- as.data.frame(idt)
colnames(idt) <- c("Treatment_Type", "Effect")
idt$se <- (lm1a$std.error)[2]
idt <- idt %>% mutate(cih90 = Effect + 1.645*se)
idt <- idt %>% mutate(cil90 = Effect - 1.645*se)
idt <- idt %>% mutate(cih95 = Effect + 1.96*se)
idt <- idt %>% mutate(cil95 = Effect - 1.96*se)

cit <- c("Whataboutism\nvs. Identity Threat", lm2a$coefficients[2]) 
cit[2] <- lapply(cit[2], as.numeric); cit <- as.data.frame(cit)
colnames(cit) <- c("Treatment_Type", "Effect")
cit$se <- (lm2a$std.error)[2]
cit <- cit %>% mutate(cih90 = Effect + 1.645*se)
cit <- cit %>% mutate(cil90 = Effect - 1.645*se)
cit <- cit %>% mutate(cih95 = Effect + 1.96*se)
cit <- cit %>% mutate(cil95 = Effect - 1.96*se)

com <- c("Whataboutism\nvs. Control", lm3a$coefficients[2]) 
com[2] <- lapply(com[2], as.numeric); com <- as.data.frame(com)
colnames(com) <- c("Treatment_Type", "Effect")
com$se <- (lm3a$std.error)[2]
com <- com %>% mutate(cih90 = Effect + 1.645*se)
com <- com %>% mutate(cil90 = Effect - 1.645*se)
com <- com %>% mutate(cih95 = Effect + 1.96*se)
com <- com %>% mutate(cil95 = Effect - 1.96*se)

#bind results by treatment type
results_age5 <- rbind(unlist(idt), unlist(cit), unlist(com))
colnames(results_age5) <- c("Treatment_Type", "Effect", "se", "cih90", "cil90", "cih95", "cil95")
results_age5 <- data.frame(results_age5)
results_age5$Effect <- as.numeric(results_age5$Effect)
results_age5$cil90 <- as.numeric(results_age5$cil90)
results_age5$cih90 <- as.numeric(results_age5$cih90)
results_age5$cil95 <- as.numeric(results_age5$cil95)
results_age5$cih95 <- as.numeric(results_age5$cih95)
results_age5$Age<- "61+"

#combine age groups and plot

age_hte <- rbind(results_age1, results_age2, results_age3, results_age4, results_age5)
age_hte$Age <- factor(age_hte$Age, levels = c("18-30", "31-40", "41-50", "51-60", "61+"))

do2 <- ggplot(age_hte, aes(x = Age, y = Effect, colour=factor(Treatment_Type, level=factor_levels))) +
  geom_point(stat="identity", position=position_dodge(width=0.5), size=3) + 
  labs(tag = "B") +
  scale_y_continuous(limits = c(-7.2, 6.6)) +
  geom_linerange(aes(ymin = cil95, ymax = cih95), 
                 position = position_dodge(width=0.5,),
                 size = 1) +
  geom_linerange(aes(ymin = cil90, ymax = cih90), 
                 position = position_dodge(width=0.5),
                 size = 2) +
  ylab("Treatment Effect") +
  xlab("Age") +
  geom_hline(yintercept=0, linetype="dashed", color="blue") +
  labs(color='Update Type') +
  theme_bw()+
  theme(
    axis.text=element_text(size=16),
    plot.title = element_text(size=25),
    axis.title=element_text(size=18),
    legend.text=element_text(size=12),
    legend.title=element_text(size=16),
    legend.key.height = unit(1.5, 'cm')
  )

####TABLE D3(C) - BLM Donation by Gender####

# women
data1_f <- subset(data1, data1$female == 1)
data2_f <- subset(data2, data2$female == 1)
data3_f <- subset(data3, data3$female == 1)

#models by treatment type
lm1a <- lm_robust(donation_amount ~ as.factor(trt) + education + income + newjob, data=data1_f, se_type="HC3")
lm2a <- lm_robust(donation_amount ~ as.factor(trt) + lr_1 + natid_covar1_1 + income + newjob, data = data2_f, se_type="HC3")
lm3a <- lm_robust(donation_amount ~ as.factor(trt) + lr_1 + unemployed, data = data3_f, se_type="HC3")

idt <- data.frame(
  Treatment_Type = "Identity Threat\nvs. Control",
  Effect = lm1a$coefficients[2],
  se = lm1a$std.error[2]
)
idt <- idt %>% 
  mutate(
    ci90h = Effect + 1.645 * se,
    ci90l = Effect - 1.645 * se,
    ci95h = Effect + 1.96 * se,
    ci95l = Effect - 1.96 * se
  )

cit <- data.frame(
  Treatment_Type = "Whataboutism\nvs. Identity Threat",
  Effect = lm2a$coefficients[2],
  se = lm2a$std.error[2]
)
cit <- cit %>% 
  mutate(
    ci90h = Effect + 1.645 * se,
    ci90l = Effect - 1.645 * se,
    ci95h = Effect + 1.96 * se,
    ci95l = Effect - 1.96 * se
  )

com <- data.frame(
  Treatment_Type = "Whataboutism\nvs. Control",
  Effect = lm3a$coefficients[2],
  se = lm3a$std.error[2]
)
com <- com %>% 
  mutate(
    ci90h = Effect + 1.645 * se,
    ci90l = Effect - 1.645 * se,
    ci95h = Effect + 1.96 * se,
    ci95l = Effect - 1.96 * se
  )

#bind results by treatment type
results_f <- rbind(idt, cit, com)
results_f$Gender <- "Female"

# men
data1_m <- subset(data1, data1$female == 0)
data2_m <- subset(data2, data2$female == 0)
data3_m <- subset(data3, data3$female == 0)

#models by treatment type
lm1a <- lm_robust(donation_amount ~ as.factor(trt) + education + income + newjob, data=data1_m, se_type="HC3")
lm2a <- lm_robust(donation_amount ~ as.factor(trt) + lr_1 + natid_covar1_1 + income + newjob, data = data2_m, se_type="HC3")
lm3a <- lm_robust(donation_amount ~ as.factor(trt) + lr_1 + unemployed, data = data3_m, se_type="HC3")

idt <- data.frame(
  Treatment_Type = "Identity Threat\nvs. Control",
  Effect = lm1a$coefficients[2],
  se = lm1a$std.error[2]
)
idt <- idt %>% 
  mutate(
    ci90h = Effect + 1.645 * se,
    ci90l = Effect - 1.645 * se,
    ci95h = Effect + 1.96 * se,
    ci95l = Effect - 1.96 * se
  )

cit <- data.frame(
  Treatment_Type = "Whataboutism\nvs. Identity Threat",
  Effect = lm2a$coefficients[2],
  se = lm2a$std.error[2]
)
cit <- cit %>% 
  mutate(
    ci90h = Effect + 1.645 * se,
    ci90l = Effect - 1.645 * se,
    ci95h = Effect + 1.96 * se,
    ci95l = Effect - 1.96 * se
  )

com <- data.frame(
  Treatment_Type = "Whataboutism\nvs. Control",
  Effect = lm3a$coefficients[2],
  se = lm3a$std.error[2]
)
com <- com %>% 
  mutate(
    ci90h = Effect + 1.645 * se,
    ci90l = Effect - 1.645 * se,
    ci95h = Effect + 1.96 * se,
    ci95l = Effect - 1.96 * se
  )

#bind results by treatment type
results_m <- rbind(idt, cit, com)
results_m$Gender <- "Male"

# Combine women/men and plot
gender_hte <- rbind(results_f, results_m)
gender_hte$Gender <- factor(gender_hte$Gender, levels = c("Female", "Male"))

do3 <- ggplot(gender_hte, aes(x = Gender, y = Effect, colour=factor(Treatment_Type, level=factor_levels))) +
  geom_point(stat="identity", position=position_dodge(width=0.3), size=3) + 
  labs(tag = "C") +
  scale_y_continuous(limits = c(-7.2, 6.6)) +
  geom_linerange(aes(ymin = ci95l, ymax = ci95h), 
                 position = position_dodge(width=0.3,),
                 size = 1) +
  geom_linerange(aes(ymin = ci90l, ymax = ci90h), 
                 position = position_dodge(width=0.3),
                 size = 2)+
  ylab("Treatment Effect") +
  xlab("Gender") +
  geom_hline(yintercept=0, linetype="dashed", color="blue") +
  labs(color='Update Type') +
  theme_bw()+
  theme(
    axis.text=element_text(size=16),
    plot.title = element_text(size=25),
    axis.title=element_text(size=18),
    legend.text=element_text(size=12),
    legend.title=element_text(size=16),
    legend.key.height = unit(1.5, 'cm')
  )


####TABLE D3(D) - BLM Donation by Pre-Treatment National Identification####

#quintile 1

data1_nid1 <- subset(data1, data1$nat_quint==1)
data2_nid1 <- subset(data2, data2$nat_quint==1)
data3_nid1 <- subset(data3, data3$nat_quint==1)

#models by treatment type
lm1a <- lm_robust(donation_amount ~ as.factor(trt) + education + income + newjob, data=data1_nid1, se_type="HC3")
lm2a <- lm_robust(donation_amount ~ as.factor(trt) + lr_1 + natid_covar1_1 + income + newjob, data = data2_nid1, se_type="HC3")
lm3a <- lm_robust(donation_amount ~ as.factor(trt) + lr_1 + unemployed, data = data3_nid1, se_type="HC3")

idt <- c("Identity Threat\nvs. Control", lm1a$coefficients[2]) 
idt[2] <- lapply(idt[2], as.numeric); idt <- as.data.frame(idt)
colnames(idt) <- c("Treatment_Type", "Effect")
idt$se <- (lm1a$std.error)[2]
idt <- idt %>% mutate(cih95 = Effect + 1.96*se)
idt <- idt %>% mutate(cil95 = Effect - 1.96*se)
idt <- idt %>% mutate(cih90 = Effect + 1.645*se)
idt <- idt %>% mutate(cil90 = Effect - 1.645*se)

cit <- c("Whataboutism\nvs. Identity Threat", lm2a$coefficients[2]) 
cit[2] <- lapply(cit[2], as.numeric); cit <- as.data.frame(cit)
colnames(cit) <- c("Treatment_Type", "Effect")
cit$se <- (lm2a$std.error)[2]
cit <- cit %>% mutate(cih95 = Effect + 1.96*se)
cit <- cit %>% mutate(cil95 = Effect - 1.96*se)
cit <- cit %>% mutate(cih90 = Effect + 1.645*se)
cit <- cit %>% mutate(cil90 = Effect - 1.645*se)

com <- c("Whataboutism\nvs. Control", lm3a$coefficients[2]) 
com[2] <- lapply(com[2], as.numeric); com <- as.data.frame(com)
colnames(com) <- c("Treatment_Type", "Effect")
com$se <- (lm3a$std.error)[2]
com <- com %>% mutate(cih95 = Effect + 1.96*se)
com <- com %>% mutate(cil95 = Effect - 1.96*se)
com <- com %>% mutate(cih90 = Effect + 1.645*se)
com <- com %>% mutate(cil90 = Effect - 1.645*se)

#bind results by treatment type
results_nid1 <- rbind(unlist(idt), unlist(cit), unlist(com))
colnames(results_nid1) <- c("Treatment_Type", "Effect", "se", "cih95", "cil95", "cih90", "cil90")
results_nid1 <- data.frame(results_nid1)
results_nid1$Effect <- as.numeric(results_nid1$Effect)
results_nid1$cil95 <- as.numeric(results_nid1$cil95)
results_nid1$cil90 <- as.numeric(results_nid1$cil90)
results_nid1$cih95 <- as.numeric(results_nid1$cih95)
results_nid1$cih90 <- as.numeric(results_nid1$cih90)
results_nid1$Quintile<- 1

#quintile 2

data1_nid2 <- subset(data1, data1$nat_quint==2)
data2_nid2 <- subset(data2, data2$nat_quint==2)
data3_nid2 <- subset(data3, data3$nat_quint==2)

#models by treatment type
lm1a <- lm_robust(donation_amount ~ as.factor(trt) + education + income + newjob, data=data1_nid2, se_type="HC3")
lm2a <- lm_robust(donation_amount ~ as.factor(trt) + lr_1 + natid_covar1_1 + income + newjob, data = data2_nid2, se_type="HC3")
lm3a <- lm_robust(donation_amount ~ as.factor(trt) + lr_1 + unemployed, data = data3_nid2, se_type="HC3")

idt <- c("Identity Threat\nvs. Control", lm1a$coefficients[2]) 
idt[2] <- lapply(idt[2], as.numeric); idt <- as.data.frame(idt)
colnames(idt) <- c("Treatment_Type", "Effect")
idt$se <- (lm1a$std.error)[2]
idt <- idt %>% mutate(cih95 = Effect + 1.96*se)
idt <- idt %>% mutate(cil95 = Effect - 1.96*se)
idt <- idt %>% mutate(cih90 = Effect + 1.645*se)
idt <- idt %>% mutate(cil90 = Effect - 1.645*se)

cit <- c("Whataboutism\nvs. Identity Threat", lm2a$coefficients[2]) 
cit[2] <- lapply(cit[2], as.numeric); cit <- as.data.frame(cit)
colnames(cit) <- c("Treatment_Type", "Effect")
cit$se <- (lm2a$std.error)[2]
cit <- cit %>% mutate(cih95 = Effect + 1.96*se)
cit <- cit %>% mutate(cil95 = Effect - 1.96*se)
cit <- cit %>% mutate(cih90 = Effect + 1.645*se)
cit <- cit %>% mutate(cil90 = Effect - 1.645*se)

com <- c("Whataboutism\nvs. Control", lm3a$coefficients[2]) 
com[2] <- lapply(com[2], as.numeric); com <- as.data.frame(com)
colnames(com) <- c("Treatment_Type", "Effect")
com$se <- (lm3a$std.error)[2]
com <- com %>% mutate(cih95 = Effect + 1.96*se)
com <- com %>% mutate(cil95 = Effect - 1.96*se)
com <- com %>% mutate(cih90 = Effect + 1.645*se)
com <- com %>% mutate(cil90 = Effect - 1.645*se)

#bind results by treatment type
results_nid2 <- rbind(unlist(idt), unlist(cit), unlist(com))
colnames(results_nid2) <- c("Treatment_Type", "Effect", "se", "cih95", "cil95", "cih90", "cil90")
results_nid2 <- data.frame(results_nid2)
results_nid2$Effect <- as.numeric(results_nid2$Effect)
results_nid2$cil95 <- as.numeric(results_nid2$cil95)
results_nid2$cil90 <- as.numeric(results_nid2$cil90)
results_nid2$cih95 <- as.numeric(results_nid2$cih95)
results_nid2$cih90 <- as.numeric(results_nid2$cih90)
results_nid2$Quintile<- 2

#quintile 3

data1_nid3 <- subset(data1, data1$nat_quint==3)
data2_nid3 <- subset(data2, data2$nat_quint==3)
data3_nid3 <- subset(data3, data3$nat_quint==3)

#models by treatment type
lm1a <- lm_robust(donation_amount ~ as.factor(trt) + education + income + newjob, data=data1_nid3, se_type="HC3")
lm2a <- lm_robust(donation_amount ~ as.factor(trt) + lr_1 + natid_covar1_1 + income + newjob, data = data2_nid3, se_type="HC3")
lm3a <- lm_robust(donation_amount ~ as.factor(trt) + lr_1 + unemployed, data = data3_nid3, se_type="HC3")

idt <- c("Identity Threat\nvs. Control", lm1a$coefficients[2]) 
idt[2] <- lapply(idt[2], as.numeric); idt <- as.data.frame(idt)
colnames(idt) <- c("Treatment_Type", "Effect")
idt$se <- (lm1a$std.error)[2]
idt <- idt %>% mutate(cih95 = Effect + 1.96*se)
idt <- idt %>% mutate(cil95 = Effect - 1.96*se)
idt <- idt %>% mutate(cih90 = Effect + 1.645*se)
idt <- idt %>% mutate(cil90 = Effect - 1.645*se)

cit <- c("Whataboutism\nvs. Identity Threat", lm2a$coefficients[2]) 
cit[2] <- lapply(cit[2], as.numeric); cit <- as.data.frame(cit)
colnames(cit) <- c("Treatment_Type", "Effect")
cit$se <- (lm2a$std.error)[2]
cit <- cit %>% mutate(cih95 = Effect + 1.96*se)
cit <- cit %>% mutate(cil95 = Effect - 1.96*se)
cit <- cit %>% mutate(cih90 = Effect + 1.645*se)
cit <- cit %>% mutate(cil90 = Effect - 1.645*se)

com <- c("Whataboutism\nvs. Control", lm3a$coefficients[2]) 
com[2] <- lapply(com[2], as.numeric); com <- as.data.frame(com)
colnames(com) <- c("Treatment_Type", "Effect")
com$se <- (lm3a$std.error)[2]
com <- com %>% mutate(cih95 = Effect + 1.96*se)
com <- com %>% mutate(cil95 = Effect - 1.96*se)
com <- com %>% mutate(cih90 = Effect + 1.645*se)
com <- com %>% mutate(cil90 = Effect - 1.645*se)

#bind results by treatment type
results_nid3 <- rbind(unlist(idt), unlist(cit), unlist(com))
colnames(results_nid3) <- c("Treatment_Type", "Effect", "se", "cih95", "cil95", "cih90", "cil90")
results_nid3 <- data.frame(results_nid3)
results_nid3$Effect <- as.numeric(results_nid3$Effect)
results_nid3$cil95 <- as.numeric(results_nid3$cil95)
results_nid3$cil90 <- as.numeric(results_nid3$cil90)
results_nid3$cih95 <- as.numeric(results_nid3$cih95)
results_nid3$cih90 <- as.numeric(results_nid3$cih90)
results_nid3$Quintile<- 3

#quintile 4

data1_nid4 <- subset(data1, data1$nat_quint==4)
data2_nid4 <- subset(data2, data2$nat_quint==4)
data3_nid4 <- subset(data3, data3$nat_quint==4)

#models by treatment type
lm1a <- lm_robust(donation_amount ~ as.factor(trt) + education + income + newjob, data=data1_nid4, se_type="HC3")
lm2a <- lm_robust(donation_amount ~ as.factor(trt) + lr_1 + natid_covar1_1 + income + newjob, data = data2_nid4, se_type="HC3")
lm3a <- lm_robust(donation_amount ~ as.factor(trt) + lr_1 + unemployed, data = data3_nid4, se_type="HC3")

idt <- c("Identity Threat\nvs. Control", lm1a$coefficients[2]) 
idt[2] <- lapply(idt[2], as.numeric); idt <- as.data.frame(idt)
colnames(idt) <- c("Treatment_Type", "Effect")
idt$se <- (lm1a$std.error)[2]
idt <- idt %>% mutate(cih95 = Effect + 1.96*se)
idt <- idt %>% mutate(cil95 = Effect - 1.96*se)
idt <- idt %>% mutate(cih90 = Effect + 1.645*se)
idt <- idt %>% mutate(cil90 = Effect - 1.645*se)

cit <- c("Whataboutism\nvs. Identity Threat", lm2a$coefficients[2]) 
cit[2] <- lapply(cit[2], as.numeric); cit <- as.data.frame(cit)
colnames(cit) <- c("Treatment_Type", "Effect")
cit$se <- (lm2a$std.error)[2]
cit <- cit %>% mutate(cih95 = Effect + 1.96*se)
cit <- cit %>% mutate(cil95 = Effect - 1.96*se)
cit <- cit %>% mutate(cih90 = Effect + 1.645*se)
cit <- cit %>% mutate(cil90 = Effect - 1.645*se)

com <- c("Whataboutism\nvs. Control", lm3a$coefficients[2]) 
com[2] <- lapply(com[2], as.numeric); com <- as.data.frame(com)
colnames(com) <- c("Treatment_Type", "Effect")
com$se <- (lm3a$std.error)[2]
com <- com %>% mutate(cih95 = Effect + 1.96*se)
com <- com %>% mutate(cil95 = Effect - 1.96*se)
com <- com %>% mutate(cih90 = Effect + 1.645*se)
com <- com %>% mutate(cil90 = Effect - 1.645*se)

#bind results by treatment type
results_nid4 <- rbind(unlist(idt), unlist(cit), unlist(com))
colnames(results_nid4) <- c("Treatment_Type", "Effect", "se", "cih95", "cil95", "cih90", "cil90")
results_nid4 <- data.frame(results_nid4)
results_nid4$Effect <- as.numeric(results_nid4$Effect)
results_nid4$cil95 <- as.numeric(results_nid4$cil95)
results_nid4$cil90 <- as.numeric(results_nid4$cil90)
results_nid4$cih95 <- as.numeric(results_nid4$cih95)
results_nid4$cih90 <- as.numeric(results_nid4$cih90)
results_nid4$Quintile<- 4

#quintile 5

data1_nid5 <- subset(data1, data1$nat_quint==5)
data2_nid5 <- subset(data2, data2$nat_quint==5)
data3_nid5 <- subset(data3, data3$nat_quint==5)

#models by treatment type
lm1a <- lm_robust(donation_amount ~ as.factor(trt) + education + income + newjob, data=data1_nid5, se_type="HC3")
lm2a <- lm_robust(donation_amount ~ as.factor(trt) + lr_1 + natid_covar1_1 + income + newjob, data = data2_nid5, se_type="HC3")
lm3a <- lm_robust(donation_amount ~ as.factor(trt) + lr_1 + unemployed, data = data3_nid5, se_type="HC3")

idt <- c("Identity Threat\nvs. Control", lm1a$coefficients[2]) 
idt[2] <- lapply(idt[2], as.numeric); idt <- as.data.frame(idt)
colnames(idt) <- c("Treatment_Type", "Effect")
idt$se <- (lm1a$std.error)[2]
idt <- idt %>% mutate(cih95 = Effect + 1.96*se)
idt <- idt %>% mutate(cil95 = Effect - 1.96*se)
idt <- idt %>% mutate(cih90 = Effect + 1.645*se)
idt <- idt %>% mutate(cil90 = Effect - 1.645*se)

cit <- c("Whataboutism\nvs. Identity Threat", lm2a$coefficients[2]) 
cit[2] <- lapply(cit[2], as.numeric); cit <- as.data.frame(cit)
colnames(cit) <- c("Treatment_Type", "Effect")
cit$se <- (lm2a$std.error)[2]
cit <- cit %>% mutate(cih95 = Effect + 1.96*se)
cit <- cit %>% mutate(cil95 = Effect - 1.96*se)
cit <- cit %>% mutate(cih90 = Effect + 1.645*se)
cit <- cit %>% mutate(cil90 = Effect - 1.645*se)

com <- c("Whataboutism\nvs. Control", lm3a$coefficients[2]) 
com[2] <- lapply(com[2], as.numeric); com <- as.data.frame(com)
colnames(com) <- c("Treatment_Type", "Effect")
com$se <- (lm3a$std.error)[2]
com <- com %>% mutate(cih95 = Effect + 1.96*se)
com <- com %>% mutate(cil95 = Effect - 1.96*se)
com <- com %>% mutate(cih90 = Effect + 1.645*se)
com <- com %>% mutate(cil90 = Effect - 1.645*se)

#bind results by treatment type
results_nid5 <- rbind(unlist(idt), unlist(cit), unlist(com))
colnames(results_nid5) <- c("Treatment_Type", "Effect", "se", "cih95", "cil95", "cih90", "cil90")
results_nid5 <- data.frame(results_nid5)
results_nid5$Effect <- as.numeric(results_nid5$Effect)
results_nid5$cil95 <- as.numeric(results_nid5$cil95)
results_nid5$cil90 <- as.numeric(results_nid5$cil90)
results_nid5$cih95 <- as.numeric(results_nid5$cih95)
results_nid5$cih90 <- as.numeric(results_nid5$cih90)
results_nid5$Quintile<- 5

#combine quintiles and plot
nid_hte <- rbind(results_nid1, results_nid2, results_nid3, results_nid4, results_nid5)
nid_hte$Effect <- as.numeric(as.character(nid_hte$Effect))

do4 <- ggplot(nid_hte, aes(x = Quintile, y = Effect, colour=factor(Treatment_Type, level=factor_levels))) +
  geom_point(stat="identity", position=position_dodge(width=0.5), size=3) + 
  labs(tag = "D") +
  scale_y_continuous(limits = c(-7.2, 6.6)) +
  geom_linerange(aes(ymin = cil95, ymax = cih95), 
                 position = position_dodge(width=0.5,),
                 size = 1) +
  geom_linerange(aes(ymin = cil90, ymax = cih90), 
                 position = position_dodge(width=0.5),
                 size = 2)+
  ylab("Treatment Effect") +
  xlab("National Identification") +
  geom_hline(yintercept=0, linetype="dashed", color="blue") +
  labs(color='Update Type') +
  theme_bw()+
  theme(
    axis.text=element_text(size=16),
    plot.title = element_text(size=25),
    axis.title=element_text(size=18),
    legend.text=element_text(size=12),
    legend.title=element_text(size=16),
    legend.key.height = unit(1.5, 'cm')
  )

####Combined Panel for Figure D3####

(do1 | do2) /
  (do3 | do4)+
  plot_layout(guides = "collect") &
  theme(legend.position = "right")

ggsave("figD3.png", width = 12, height = 10, dpi = 300)

####TABLE D4 - HTEs for National Identification####
####TABLE D4(A) - National Identification by Left/Right Self-Placement####

#Left-wing
data1_left <- subset(data1, data1$lr_1 < 5)
data2_left <- subset(data2, data2$lr_1 < 5)
data3_left <- subset(data3, data3$lr_1 < 5)

#models by treatment type
lm1a <- lm_robust(natid_pca ~ as.factor(trt) + education + income + newjob, data=data1_left, se_type="HC3")
lm2a <- lm_robust(natid_pca ~ as.factor(trt) + lr_1 + natid_covar1_1 + income + newjob, data = data2_left, se_type="HC3")
lm3a <- lm_robust(natid_pca ~ as.factor(trt) + lr_1 + unemployed, data = data3_left, se_type="HC3")

idt <- c("Identity Threat\nvs. Control", lm1a$coefficients[2])
idt[2] <- lapply(idt[2], as.numeric); idt <- as.data.frame(idt)
colnames(idt) <- c("Treatment_Type", "Effect")
idt$se <- lm1a$std.error[2]
idt <- idt %>% 
  mutate(
    cih_95 = Effect + z_95 * se,
    cil_95 = Effect - z_95 * se,
    cih_90 = Effect + z_90 * se,
    cil_90 = Effect - z_90 * se
  )

cit <- c("Whataboutism\nvs. Identity Threat", lm2a$coefficients[2])
cit[2] <- lapply(cit[2], as.numeric); cit <- as.data.frame(cit)
colnames(cit) <- c("Treatment_Type", "Effect")
cit$se <- lm2a$std.error[2]
cit <- cit %>% 
  mutate(
    cih_95 = Effect + z_95 * se,
    cil_95 = Effect - z_95 * se,
    cih_90 = Effect + z_90 * se,
    cil_90 = Effect - z_90 * se
  )

com <- c("Whataboutism\nvs. Control", lm3a$coefficients[2])
com[2] <- lapply(com[2], as.numeric); com <- as.data.frame(com)
colnames(com) <- c("Treatment_Type", "Effect")
com$se <- lm3a$std.error[2]
com <- com %>% 
  mutate(
    cih_95 = Effect + z_95 * se,
    cil_95 = Effect - z_95 * se,
    cih_90 = Effect + z_90 * se,
    cil_90 = Effect - z_90 * se
  )

#bind results by treatment type
results_left <- rbind(unlist(idt), unlist(cit), unlist(com))
colnames(results_left) <- c("Treatment_Type", "Effect", "se", "cih_95", "cil_95", "cih_90", "cil_90")
results_left <- data.frame(results_left)
results_left$Effect <- as.numeric(results_left$Effect)
results_left$cil_95 <- as.numeric(results_left$cil_95)
results_left$cih_95 <- as.numeric(results_left$cih_95)
results_left$cil_90 <- as.numeric(results_left$cil_90)
results_left$cih_90 <- as.numeric(results_left$cih_90)
results_left$HTE <- "Left"

#Right-wing
data1_right <- subset(data1, data1$lr_1 >= 6)
data2_right <- subset(data2, data2$lr_1 >= 6)
data3_right <- subset(data3, data3$lr_1 >= 6)

#models by treatment type
lm1a <- lm_robust(natid_pca ~ as.factor(trt) + education + income + newjob, data=data1_right, se_type="HC3")
lm2a <- lm_robust(natid_pca ~ as.factor(trt) + lr_1 + natid_covar1_1 + income + newjob, data = data2_right, se_type="HC3")
lm3a <- lm_robust(natid_pca ~ as.factor(trt) + lr_1 + unemployed, data = data3_right, se_type="HC3")

idt <- c("Identity Threat\nvs. Control", lm1a$coefficients[2])
idt[2] <- lapply(idt[2], as.numeric); idt <- as.data.frame(idt)
colnames(idt) <- c("Treatment_Type", "Effect")
idt$se <- lm1a$std.error[2]
idt <- idt %>% 
  mutate(
    cih_95 = Effect + z_95 * se,
    cil_95 = Effect - z_95 * se,
    cih_90 = Effect + z_90 * se,
    cil_90 = Effect - z_90 * se
  )

cit <- c("Whataboutism\nvs. Identity Threat", lm2a$coefficients[2])
cit[2] <- lapply(cit[2], as.numeric); cit <- as.data.frame(cit)
colnames(cit) <- c("Treatment_Type", "Effect")
cit$se <- lm2a$std.error[2]
cit <- cit %>% 
  mutate(
    cih_95 = Effect + z_95 * se,
    cil_95 = Effect - z_95 * se,
    cih_90 = Effect + z_90 * se,
    cil_90 = Effect - z_90 * se
  )

com <- c("Whataboutism\nvs. Control", lm3a$coefficients[2])
com[2] <- lapply(com[2], as.numeric); com <- as.data.frame(com)
colnames(com) <- c("Treatment_Type", "Effect")
com$se <- lm3a$std.error[2]
com <- com %>% 
  mutate(
    cih_95 = Effect + z_95 * se,
    cil_95 = Effect - z_95 * se,
    cih_90 = Effect + z_90 * se,
    cil_90 = Effect - z_90 * se
  )

#bind results by treatment type
results_right <- rbind(unlist(idt), unlist(cit), unlist(com))
colnames(results_right) <- c("Treatment_Type", "Effect", "se", "cih_95", "cil_95", "cih_90", "cil_90")
results_right <- data.frame(results_right)
results_right$Effect <- as.numeric(results_right$Effect)
results_right$cil_95 <- as.numeric(results_right$cil_95)
results_right$cih_95 <- as.numeric(results_right$cih_95)
results_right$cil_90 <- as.numeric(results_right$cil_90)
results_right$cih_90 <- as.numeric(results_right$cih_90)
results_right$HTE <- "Right"

#bind left and right results and plot
lr_hte <- rbind(results_left, results_right)
lr_hte$HTE <- factor(lr_hte$HTE, levels = c("Left", "Right"))

ni1 <- ggplot(lr_hte, aes(x = HTE, y = Effect, colour=factor(Treatment_Type, level=factor_levels))) +
  geom_point(stat="identity", position=position_dodge(width=0.3), size=3) + 
  labs(tag = "A") +
  scale_y_continuous(limits = c(-0.55, 0.35)) +
  geom_linerange(aes(ymin = cil_95, ymax = cih_95), 
                 position = position_dodge(width=0.3,),
                 size = 1) +
  geom_linerange(aes(ymin = cil_90, ymax = cih_90), 
                 position = position_dodge(width=0.3),
                 size = 2)+
  ylab("Treatment Effect") +
  xlab("Left-Right Self-Placement") +
  geom_hline(yintercept=0, linetype="dashed", color="blue") +
  labs(color='Update Type') +
  theme_bw()+
  theme(
    axis.text=element_text(size=16),
    plot.title = element_text(size=25),
    axis.title=element_text(size=18),
    legend.text=element_text(size=12),
    legend.title=element_text(size=16),
    legend.key.height = unit(1.5, 'cm')
  )

####TABLE D4(B) - National Identification by Age####

#age group 1

data1_age1 <- subset(data1, data1$age<31)
data2_age1 <- subset(data2, data2$age<31)
data3_age1 <- subset(data3, data3$age<31)

#models by treatment type
lm1a <- lm_robust(natid_pca ~ as.factor(trt) + education + income + newjob, data=data1_age1, se_type="HC3")
lm2a <- lm_robust(natid_pca ~ as.factor(trt) + lr_1  + natid_covar1_1 + income + newjob, data=data2_age1, se_type="HC3")
lm3a <- lm_robust(natid_pca ~ as.factor(trt) + lr_1 + unemployed, data=data3_age1, se_type="HC3")

idt <- c("Identity Threat\nvs. Control", lm1a$coefficients[2]) 
idt[2] <- lapply(idt[2], as.numeric); idt <- as.data.frame(idt)
colnames(idt) <- c("Treatment_Type", "Effect")
idt$se <- (lm1a$std.error)[2]
idt <- idt %>% mutate(cih90 = Effect + 1.645*se)
idt <- idt %>% mutate(cil90 = Effect - 1.645*se)
idt <- idt %>% mutate(cih95 = Effect + 1.96*se)
idt <- idt %>% mutate(cil95 = Effect - 1.96*se)

cit <- c("Whataboutism\nvs. Identity Threat", lm2a$coefficients[2]) 
cit[2] <- lapply(cit[2], as.numeric); cit <- as.data.frame(cit)
colnames(cit) <- c("Treatment_Type", "Effect")
cit$se <- (lm2a$std.error)[2]
cit <- cit %>% mutate(cih90 = Effect + 1.645*se)
cit <- cit %>% mutate(cil90 = Effect - 1.645*se)
cit <- cit %>% mutate(cih95 = Effect + 1.96*se)
cit <- cit %>% mutate(cil95 = Effect - 1.96*se)

com <- c("Whataboutism\nvs. Control", lm3a$coefficients[2]) 
com[2] <- lapply(com[2], as.numeric); com <- as.data.frame(com)
colnames(com) <- c("Treatment_Type", "Effect")
com$se <- (lm3a$std.error)[2]
com <- com %>% mutate(cih90 = Effect + 1.645*se)
com <- com %>% mutate(cil90 = Effect - 1.645*se)
com <- com %>% mutate(cih95 = Effect + 1.96*se)
com <- com %>% mutate(cil95 = Effect - 1.96*se)

#bind results by treatment type
results_age1 <- rbind(unlist(idt), unlist(cit), unlist(com))
colnames(results_age1) <- c("Treatment_Type", "Effect", "se", "cih90", "cil90", "cih95", "cil95")
results_age1 <- data.frame(results_age1)
results_age1$Effect <- as.numeric(results_age1$Effect)
results_age1$cil90 <- as.numeric(results_age1$cil90)
results_age1$cih90 <- as.numeric(results_age1$cih90)
results_age1$cil95 <- as.numeric(results_age1$cil95)
results_age1$cih95 <- as.numeric(results_age1$cih95)

results_age1$Age<- "18-30"

#age group 2

data1_age2 <- subset(data1, data1$age>30 & data1$age<41)
data2_age2 <- subset(data2, data2$age>30 & data2$age<41)
data3_age2 <- subset(data3, data3$age>30 & data3$age<41)

#models by treatment type
lm1a <- lm_robust(natid_pca ~ as.factor(trt) + education + income + newjob, data=data1_age2, se_type="HC3")
lm2a <- lm_robust(natid_pca ~ as.factor(trt) + lr_1  + natid_covar1_1 + income + newjob, data=data2_age2, se_type="HC3")
lm3a <- lm_robust(natid_pca ~ as.factor(trt) + lr_1 + unemployed, data=data3_age2, se_type="HC3")

idt <- c("Identity Threat\nvs. Control", lm1a$coefficients[2]) 
idt[2] <- lapply(idt[2], as.numeric); idt <- as.data.frame(idt)
colnames(idt) <- c("Treatment_Type", "Effect")
idt$se <- (lm1a$std.error)[2]
idt <- idt %>% mutate(cih90 = Effect + 1.645*se)
idt <- idt %>% mutate(cil90 = Effect - 1.645*se)
idt <- idt %>% mutate(cih95 = Effect + 1.96*se)
idt <- idt %>% mutate(cil95 = Effect - 1.96*se)

cit <- c("Whataboutism\nvs. Identity Threat", lm2a$coefficients[2]) 
cit[2] <- lapply(cit[2], as.numeric); cit <- as.data.frame(cit)
colnames(cit) <- c("Treatment_Type", "Effect")
cit$se <- (lm2a$std.error)[2]
cit <- cit %>% mutate(cih90 = Effect + 1.645*se)
cit <- cit %>% mutate(cil90 = Effect - 1.645*se)
cit <- cit %>% mutate(cih95 = Effect + 1.96*se)
cit <- cit %>% mutate(cil95 = Effect - 1.96*se)

com <- c("Whataboutism\nvs. Control", lm3a$coefficients[2]) 
com[2] <- lapply(com[2], as.numeric); com <- as.data.frame(com)
colnames(com) <- c("Treatment_Type", "Effect")
com$se <- (lm3a$std.error)[2]
com <- com %>% mutate(cih90 = Effect + 1.645*se)
com <- com %>% mutate(cil90 = Effect - 1.645*se)
com <- com %>% mutate(cih95 = Effect + 1.96*se)
com <- com %>% mutate(cil95 = Effect - 1.96*se)

#bind results by treatment type
results_age2 <- rbind(unlist(idt), unlist(cit), unlist(com))
colnames(results_age2) <- c("Treatment_Type", "Effect", "se", "cih90", "cil90", "cih95", "cil95")
results_age2 <- data.frame(results_age2)
results_age2$Effect <- as.numeric(results_age2$Effect)
results_age2$cil90 <- as.numeric(results_age2$cil90)
results_age2$cih90 <- as.numeric(results_age2$cih90)
results_age2$cil95 <- as.numeric(results_age2$cil95)
results_age2$cih95 <- as.numeric(results_age2$cih95)
results_age2$Age<- "31-40"

#age group 3

data1_age3 <- subset(data1, data1$age>40 & data1$age<51)
data2_age3 <- subset(data2, data2$age>40 & data2$age<51)
data3_age3 <- subset(data3, data3$age>40 & data3$age<51)

#models by treatment type
lm1a <- lm_robust(natid_pca ~ as.factor(trt) + education + income + newjob, data=data1_age3, se_type="HC3")
lm2a <- lm_robust(natid_pca ~ as.factor(trt) + lr_1  + natid_covar1_1 + income + newjob, data=data2_age3, se_type="HC3")
lm3a <- lm_robust(natid_pca ~ as.factor(trt) + lr_1 + unemployed, data=data3_age3, se_type="HC3")

idt <- c("Identity Threat\nvs. Control", lm1a$coefficients[2]) 
idt[2] <- lapply(idt[2], as.numeric); idt <- as.data.frame(idt)
colnames(idt) <- c("Treatment_Type", "Effect")
idt$se <- (lm1a$std.error)[2]
idt <- idt %>% mutate(cih90 = Effect + 1.645*se)
idt <- idt %>% mutate(cil90 = Effect - 1.645*se)
idt <- idt %>% mutate(cih95 = Effect + 1.96*se)
idt <- idt %>% mutate(cil95 = Effect - 1.96*se)

cit <- c("Whataboutism\nvs. Identity Threat", lm2a$coefficients[2]) 
cit[2] <- lapply(cit[2], as.numeric); cit <- as.data.frame(cit)
colnames(cit) <- c("Treatment_Type", "Effect")
cit$se <- (lm2a$std.error)[2]
cit <- cit %>% mutate(cih90 = Effect + 1.645*se)
cit <- cit %>% mutate(cil90 = Effect - 1.645*se)
cit <- cit %>% mutate(cih95 = Effect + 1.96*se)
cit <- cit %>% mutate(cil95 = Effect - 1.96*se)

com <- c("Whataboutism\nvs. Control", lm3a$coefficients[2]) 
com[2] <- lapply(com[2], as.numeric); com <- as.data.frame(com)
colnames(com) <- c("Treatment_Type", "Effect")
com$se <- (lm3a$std.error)[2]
com <- com %>% mutate(cih90 = Effect + 1.645*se)
com <- com %>% mutate(cil90 = Effect - 1.645*se)
com <- com %>% mutate(cih95 = Effect + 1.96*se)
com <- com %>% mutate(cil95 = Effect - 1.96*se)

#bind results by treatment type
results_age3 <- rbind(unlist(idt), unlist(cit), unlist(com))
colnames(results_age3) <- c("Treatment_Type", "Effect", "se", "cih90", "cil90", "cih95", "cil95")
results_age3 <- data.frame(results_age3)
results_age3$Effect <- as.numeric(results_age3$Effect)
results_age3$cil90 <- as.numeric(results_age3$cil90)
results_age3$cih90 <- as.numeric(results_age3$cih90)
results_age3$cil95 <- as.numeric(results_age3$cil95)
results_age3$cih95 <- as.numeric(results_age3$cih95)
results_age3$Age<- "41-50"

#age group 4

data1_age4 <- subset(data1, data1$age>50 & data1$age<61)
data2_age4 <- subset(data2, data2$age>50 & data2$age<61)
data3_age4 <- subset(data3, data3$age>50 & data3$age<61)

#models by treatment type
lm1a <- lm_robust(natid_pca ~ as.factor(trt) + education + income + newjob, data=data1_age4, se_type="HC3")
lm2a <- lm_robust(natid_pca ~ as.factor(trt) + lr_1  + natid_covar1_1 + income + newjob, data=data2_age4, se_type="HC3")
lm3a <- lm_robust(natid_pca ~ as.factor(trt) + lr_1 + unemployed, data=data3_age4, se_type="HC3")

idt <- c("Identity Threat\nvs. Control", lm1a$coefficients[2]) 
idt[2] <- lapply(idt[2], as.numeric); idt <- as.data.frame(idt)
colnames(idt) <- c("Treatment_Type", "Effect")
idt$se <- (lm1a$std.error)[2]
idt <- idt %>% mutate(cih90 = Effect + 1.645*se)
idt <- idt %>% mutate(cil90 = Effect - 1.645*se)
idt <- idt %>% mutate(cih95 = Effect + 1.96*se)
idt <- idt %>% mutate(cil95 = Effect - 1.96*se)

cit <- c("Whataboutism\nvs. Identity Threat", lm2a$coefficients[2]) 
cit[2] <- lapply(cit[2], as.numeric); cit <- as.data.frame(cit)
colnames(cit) <- c("Treatment_Type", "Effect")
cit$se <- (lm2a$std.error)[2]
cit <- cit %>% mutate(cih90 = Effect + 1.645*se)
cit <- cit %>% mutate(cil90 = Effect - 1.645*se)
cit <- cit %>% mutate(cih95 = Effect + 1.96*se)
cit <- cit %>% mutate(cil95 = Effect - 1.96*se)

com <- c("Whataboutism\nvs. Control", lm3a$coefficients[2]) 
com[2] <- lapply(com[2], as.numeric); com <- as.data.frame(com)
colnames(com) <- c("Treatment_Type", "Effect")
com$se <- (lm3a$std.error)[2]
com <- com %>% mutate(cih90 = Effect + 1.645*se)
com <- com %>% mutate(cil90 = Effect - 1.645*se)
com <- com %>% mutate(cih95 = Effect + 1.96*se)
com <- com %>% mutate(cil95 = Effect - 1.96*se)

#bind results by treatment type
results_age4 <- rbind(unlist(idt), unlist(cit), unlist(com))
colnames(results_age4) <- c("Treatment_Type", "Effect", "se", "cih90", "cil90", "cih95", "cil95")
results_age4 <- data.frame(results_age4)
results_age4$Effect <- as.numeric(results_age4$Effect)
results_age4$cil90 <- as.numeric(results_age4$cil90)
results_age4$cih90 <- as.numeric(results_age4$cih90)
results_age4$cil95 <- as.numeric(results_age4$cil95)
results_age4$cih95 <- as.numeric(results_age4$cih95)
results_age4$Age<- "51-60"

#age group 5

data1_age5 <- subset(data1, data1$age>60)
data2_age5 <- subset(data2, data2$age>60)
data3_age5 <- subset(data3, data3$age>60)

#models by treatment type
lm1a <- lm_robust(natid_pca ~ as.factor(trt) + education + income + newjob, data=data1_age5, se_type="HC3")
lm2a <- lm_robust(natid_pca ~ as.factor(trt) + lr_1  + natid_covar1_1 + income + newjob, data=data2_age5, se_type="HC3")
lm3a <- lm_robust(natid_pca ~ as.factor(trt) + lr_1 + unemployed, data=data3_age5, se_type="HC3")

idt <- c("Identity Threat\nvs. Control", lm1a$coefficients[2]) 
idt[2] <- lapply(idt[2], as.numeric); idt <- as.data.frame(idt)
colnames(idt) <- c("Treatment_Type", "Effect")
idt$se <- (lm1a$std.error)[2]
idt <- idt %>% mutate(cih90 = Effect + 1.645*se)
idt <- idt %>% mutate(cil90 = Effect - 1.645*se)
idt <- idt %>% mutate(cih95 = Effect + 1.96*se)
idt <- idt %>% mutate(cil95 = Effect - 1.96*se)

cit <- c("Whataboutism\nvs. Identity Threat", lm2a$coefficients[2]) 
cit[2] <- lapply(cit[2], as.numeric); cit <- as.data.frame(cit)
colnames(cit) <- c("Treatment_Type", "Effect")
cit$se <- (lm2a$std.error)[2]
cit <- cit %>% mutate(cih90 = Effect + 1.645*se)
cit <- cit %>% mutate(cil90 = Effect - 1.645*se)
cit <- cit %>% mutate(cih95 = Effect + 1.96*se)
cit <- cit %>% mutate(cil95 = Effect - 1.96*se)

com <- c("Whataboutism\nvs. Control", lm3a$coefficients[2]) 
com[2] <- lapply(com[2], as.numeric); com <- as.data.frame(com)
colnames(com) <- c("Treatment_Type", "Effect")
com$se <- (lm3a$std.error)[2]
com <- com %>% mutate(cih90 = Effect + 1.645*se)
com <- com %>% mutate(cil90 = Effect - 1.645*se)
com <- com %>% mutate(cih95 = Effect + 1.96*se)
com <- com %>% mutate(cil95 = Effect - 1.96*se)

#bind results by treatment type
results_age5 <- rbind(unlist(idt), unlist(cit), unlist(com))
colnames(results_age5) <- c("Treatment_Type", "Effect", "se", "cih90", "cil90", "cih95", "cil95")
results_age5 <- data.frame(results_age5)
results_age5$Effect <- as.numeric(results_age5$Effect)
results_age5$cil90 <- as.numeric(results_age5$cil90)
results_age5$cih90 <- as.numeric(results_age5$cih90)
results_age5$cil95 <- as.numeric(results_age5$cil95)
results_age5$cih95 <- as.numeric(results_age5$cih95)
results_age5$Age<- "61+"

#combine age groups and plot

age_hte <- rbind(results_age1, results_age2, results_age3, results_age4, results_age5)
age_hte$Age <- factor(age_hte$Age, levels = c("18-30", "31-40", "41-50", "51-60", "61+"))

ni2 <- ggplot(age_hte, aes(x = Age, y = Effect, colour=factor(Treatment_Type, level=factor_levels))) +
  geom_point(stat="identity", position=position_dodge(width=0.5), size=3) + 
  labs(tag = "B") +
  scale_y_continuous(limits = c(-0.55, 0.35)) +
  geom_linerange(aes(ymin = cil95, ymax = cih95), 
                 position = position_dodge(width=0.5,),
                 size = 1) +
  geom_linerange(aes(ymin = cil90, ymax = cih90), 
                 position = position_dodge(width=0.5),
                 size = 2) +
  ylab("Treatment Effect") +
  xlab("Age") +
  geom_hline(yintercept=0, linetype="dashed", color="blue") +
  labs(color='Update Type') +
  theme_bw()+
  theme(
    axis.text=element_text(size=16),
    plot.title = element_text(size=25),
    axis.title=element_text(size=18),
    legend.text=element_text(size=12),
    legend.title=element_text(size=16),
    legend.key.height = unit(1.5, 'cm')
  )

####TABLE D4(C) - National Identification by Gender####

# women
data1_f <- subset(data1, data1$female == 1)
data2_f <- subset(data2, data2$female == 1)
data3_f <- subset(data3, data3$female == 1)

#models by treatment type
lm1a <- lm_robust(natid_pca ~ as.factor(trt) + education + income + newjob, data=data1_f, se_type="HC3")
lm2a <- lm_robust(natid_pca ~ as.factor(trt) + lr_1 + natid_covar1_1 + income + newjob, data = data2_f, se_type="HC3")
lm3a <- lm_robust(natid_pca ~ as.factor(trt) + lr_1 + unemployed, data = data3_f, se_type="HC3")

idt <- data.frame(
  Treatment_Type = "Identity Threat\nvs. Control",
  Effect = lm1a$coefficients[2],
  se = lm1a$std.error[2]
)
idt <- idt %>% 
  mutate(
    ci90h = Effect + 1.645 * se,
    ci90l = Effect - 1.645 * se,
    ci95h = Effect + 1.96 * se,
    ci95l = Effect - 1.96 * se
  )

cit <- data.frame(
  Treatment_Type = "Whataboutism\nvs. Identity Threat",
  Effect = lm2a$coefficients[2],
  se = lm2a$std.error[2]
)
cit <- cit %>% 
  mutate(
    ci90h = Effect + 1.645 * se,
    ci90l = Effect - 1.645 * se,
    ci95h = Effect + 1.96 * se,
    ci95l = Effect - 1.96 * se
  )

com <- data.frame(
  Treatment_Type = "Whataboutism\nvs. Control",
  Effect = lm3a$coefficients[2],
  se = lm3a$std.error[2]
)
com <- com %>% 
  mutate(
    ci90h = Effect + 1.645 * se,
    ci90l = Effect - 1.645 * se,
    ci95h = Effect + 1.96 * se,
    ci95l = Effect - 1.96 * se
  )

#bind results by treatment type
results_f <- rbind(idt, cit, com)
results_f$Gender <- "Female"

# men
data1_m <- subset(data1, data1$female == 0)
data2_m <- subset(data2, data2$female == 0)
data3_m <- subset(data3, data3$female == 0)

#models by treatment type
lm1a <- lm_robust(natid_pca ~ as.factor(trt) + education + income + newjob, data=data1_m, se_type="HC3")
lm2a <- lm_robust(natid_pca ~ as.factor(trt) + lr_1 + natid_covar1_1 + income + newjob, data = data2_m, se_type="HC3")
lm3a <- lm_robust(natid_pca ~ as.factor(trt) + lr_1 + unemployed, data = data3_m, se_type="HC3")

idt <- data.frame(
  Treatment_Type = "Identity Threat\nvs. Control",
  Effect = lm1a$coefficients[2],
  se = lm1a$std.error[2]
)
idt <- idt %>% 
  mutate(
    ci90h = Effect + 1.645 * se,
    ci90l = Effect - 1.645 * se,
    ci95h = Effect + 1.96 * se,
    ci95l = Effect - 1.96 * se
  )

cit <- data.frame(
  Treatment_Type = "Whataboutism\nvs. Identity Threat",
  Effect = lm2a$coefficients[2],
  se = lm2a$std.error[2]
)
cit <- cit %>% 
  mutate(
    ci90h = Effect + 1.645 * se,
    ci90l = Effect - 1.645 * se,
    ci95h = Effect + 1.96 * se,
    ci95l = Effect - 1.96 * se
  )

com <- data.frame(
  Treatment_Type = "Whataboutism\nvs. Control",
  Effect = lm3a$coefficients[2],
  se = lm3a$std.error[2]
)
com <- com %>% 
  mutate(
    ci90h = Effect + 1.645 * se,
    ci90l = Effect - 1.645 * se,
    ci95h = Effect + 1.96 * se,
    ci95l = Effect - 1.96 * se
  )

#bind results by treatment type
results_m <- rbind(idt, cit, com)
results_m$Gender <- "Male"

# Combine women/men and plot
gender_hte <- rbind(results_f, results_m)
gender_hte$Gender <- factor(gender_hte$Gender, levels = c("Female", "Male"))

ni3 <- ggplot(gender_hte, aes(x = Gender, y = Effect, colour=factor(Treatment_Type, level=factor_levels))) +
  geom_point(stat="identity", position=position_dodge(width=0.3), size=3) + 
  labs(tag = "C") +
  scale_y_continuous(limits = c(-0.55, 0.35)) +
  geom_linerange(aes(ymin = ci95l, ymax = ci95h), 
                 position = position_dodge(width=0.3,),
                 size = 1) +
  geom_linerange(aes(ymin = ci90l, ymax = ci90h), 
                 position = position_dodge(width=0.3),
                 size = 2)+
  ylab("Treatment Effect") +
  xlab("Gender") +
  geom_hline(yintercept=0, linetype="dashed", color="blue") +
  labs(color='Update Type') +
  theme_bw()+
  theme(
    axis.text=element_text(size=16),
    plot.title = element_text(size=25),
    axis.title=element_text(size=18),
    legend.text=element_text(size=12),
    legend.title=element_text(size=16),
    legend.key.height = unit(1.5, 'cm')
  )


####TABLE D4(D) - National Identification by Pre-Treatment National Identification####

#quintile 1

data1_nid1 <- subset(data1, data1$nat_quint==1)
data2_nid1 <- subset(data2, data2$nat_quint==1)
data3_nid1 <- subset(data3, data3$nat_quint==1)

#models by treatment type
lm1a <- lm_robust(natid_pca ~ as.factor(trt) + education + income + newjob, data=data1_nid1, se_type="HC3")
lm2a <- lm_robust(natid_pca ~ as.factor(trt) + lr_1 + natid_covar1_1 + income + newjob, data = data2_nid1, se_type="HC3")
lm3a <- lm_robust(natid_pca ~ as.factor(trt) + lr_1 + unemployed, data = data3_nid1, se_type="HC3")

idt <- c("Identity Threat\nvs. Control", lm1a$coefficients[2]) 
idt[2] <- lapply(idt[2], as.numeric); idt <- as.data.frame(idt)
colnames(idt) <- c("Treatment_Type", "Effect")
idt$se <- (lm1a$std.error)[2]
idt <- idt %>% mutate(cih95 = Effect + 1.96*se)
idt <- idt %>% mutate(cil95 = Effect - 1.96*se)
idt <- idt %>% mutate(cih90 = Effect + 1.645*se)
idt <- idt %>% mutate(cil90 = Effect - 1.645*se)

cit <- c("Whataboutism\nvs. Identity Threat", lm2a$coefficients[2]) 
cit[2] <- lapply(cit[2], as.numeric); cit <- as.data.frame(cit)
colnames(cit) <- c("Treatment_Type", "Effect")
cit$se <- (lm2a$std.error)[2]
cit <- cit %>% mutate(cih95 = Effect + 1.96*se)
cit <- cit %>% mutate(cil95 = Effect - 1.96*se)
cit <- cit %>% mutate(cih90 = Effect + 1.645*se)
cit <- cit %>% mutate(cil90 = Effect - 1.645*se)

com <- c("Whataboutism\nvs. Control", lm3a$coefficients[2]) 
com[2] <- lapply(com[2], as.numeric); com <- as.data.frame(com)
colnames(com) <- c("Treatment_Type", "Effect")
com$se <- (lm3a$std.error)[2]
com <- com %>% mutate(cih95 = Effect + 1.96*se)
com <- com %>% mutate(cil95 = Effect - 1.96*se)
com <- com %>% mutate(cih90 = Effect + 1.645*se)
com <- com %>% mutate(cil90 = Effect - 1.645*se)

#bind results by treatment type
results_nid1 <- rbind(unlist(idt), unlist(cit), unlist(com))
colnames(results_nid1) <- c("Treatment_Type", "Effect", "se", "cih95", "cil95", "cih90", "cil90")
results_nid1 <- data.frame(results_nid1)
results_nid1$Effect <- as.numeric(results_nid1$Effect)
results_nid1$cil95 <- as.numeric(results_nid1$cil95)
results_nid1$cil90 <- as.numeric(results_nid1$cil90)
results_nid1$cih95 <- as.numeric(results_nid1$cih95)
results_nid1$cih90 <- as.numeric(results_nid1$cih90)
results_nid1$Quintile<- 1

#quintile 2

data1_nid2 <- subset(data1, data1$nat_quint==2)
data2_nid2 <- subset(data2, data2$nat_quint==2)
data3_nid2 <- subset(data3, data3$nat_quint==2)

#models by treatment type
lm1a <- lm_robust(natid_pca ~ as.factor(trt) + education + income + newjob, data=data1_nid2, se_type="HC3")
lm2a <- lm_robust(natid_pca ~ as.factor(trt) + lr_1 + natid_covar1_1 + income + newjob, data = data2_nid2, se_type="HC3")
lm3a <- lm_robust(natid_pca ~ as.factor(trt) + lr_1 + unemployed, data = data3_nid2, se_type="HC3")

idt <- c("Identity Threat\nvs. Control", lm1a$coefficients[2]) 
idt[2] <- lapply(idt[2], as.numeric); idt <- as.data.frame(idt)
colnames(idt) <- c("Treatment_Type", "Effect")
idt$se <- (lm1a$std.error)[2]
idt <- idt %>% mutate(cih95 = Effect + 1.96*se)
idt <- idt %>% mutate(cil95 = Effect - 1.96*se)
idt <- idt %>% mutate(cih90 = Effect + 1.645*se)
idt <- idt %>% mutate(cil90 = Effect - 1.645*se)

cit <- c("Whataboutism\nvs. Identity Threat", lm2a$coefficients[2]) 
cit[2] <- lapply(cit[2], as.numeric); cit <- as.data.frame(cit)
colnames(cit) <- c("Treatment_Type", "Effect")
cit$se <- (lm2a$std.error)[2]
cit <- cit %>% mutate(cih95 = Effect + 1.96*se)
cit <- cit %>% mutate(cil95 = Effect - 1.96*se)
cit <- cit %>% mutate(cih90 = Effect + 1.645*se)
cit <- cit %>% mutate(cil90 = Effect - 1.645*se)

com <- c("Whataboutism\nvs. Control", lm3a$coefficients[2]) 
com[2] <- lapply(com[2], as.numeric); com <- as.data.frame(com)
colnames(com) <- c("Treatment_Type", "Effect")
com$se <- (lm3a$std.error)[2]
com <- com %>% mutate(cih95 = Effect + 1.96*se)
com <- com %>% mutate(cil95 = Effect - 1.96*se)
com <- com %>% mutate(cih90 = Effect + 1.645*se)
com <- com %>% mutate(cil90 = Effect - 1.645*se)

#bind results by treatment type
results_nid2 <- rbind(unlist(idt), unlist(cit), unlist(com))
colnames(results_nid2) <- c("Treatment_Type", "Effect", "se", "cih95", "cil95", "cih90", "cil90")
results_nid2 <- data.frame(results_nid2)
results_nid2$Effect <- as.numeric(results_nid2$Effect)
results_nid2$cil95 <- as.numeric(results_nid2$cil95)
results_nid2$cil90 <- as.numeric(results_nid2$cil90)
results_nid2$cih95 <- as.numeric(results_nid2$cih95)
results_nid2$cih90 <- as.numeric(results_nid2$cih90)
results_nid2$Quintile<- 2

#quintile 3

data1_nid3 <- subset(data1, data1$nat_quint==3)
data2_nid3 <- subset(data2, data2$nat_quint==3)
data3_nid3 <- subset(data3, data3$nat_quint==3)

#models by treatment type
lm1a <- lm_robust(natid_pca ~ as.factor(trt) + education + income + newjob, data=data1_nid3, se_type="HC3")
lm2a <- lm_robust(natid_pca ~ as.factor(trt) + lr_1 + natid_covar1_1 + income + newjob, data = data2_nid3, se_type="HC3")
lm3a <- lm_robust(natid_pca ~ as.factor(trt) + lr_1 + unemployed, data = data3_nid3, se_type="HC3")

idt <- c("Identity Threat\nvs. Control", lm1a$coefficients[2]) 
idt[2] <- lapply(idt[2], as.numeric); idt <- as.data.frame(idt)
colnames(idt) <- c("Treatment_Type", "Effect")
idt$se <- (lm1a$std.error)[2]
idt <- idt %>% mutate(cih95 = Effect + 1.96*se)
idt <- idt %>% mutate(cil95 = Effect - 1.96*se)
idt <- idt %>% mutate(cih90 = Effect + 1.645*se)
idt <- idt %>% mutate(cil90 = Effect - 1.645*se)

cit <- c("Whataboutism\nvs. Identity Threat", lm2a$coefficients[2]) 
cit[2] <- lapply(cit[2], as.numeric); cit <- as.data.frame(cit)
colnames(cit) <- c("Treatment_Type", "Effect")
cit$se <- (lm2a$std.error)[2]
cit <- cit %>% mutate(cih95 = Effect + 1.96*se)
cit <- cit %>% mutate(cil95 = Effect - 1.96*se)
cit <- cit %>% mutate(cih90 = Effect + 1.645*se)
cit <- cit %>% mutate(cil90 = Effect - 1.645*se)

com <- c("Whataboutism\nvs. Control", lm3a$coefficients[2]) 
com[2] <- lapply(com[2], as.numeric); com <- as.data.frame(com)
colnames(com) <- c("Treatment_Type", "Effect")
com$se <- (lm3a$std.error)[2]
com <- com %>% mutate(cih95 = Effect + 1.96*se)
com <- com %>% mutate(cil95 = Effect - 1.96*se)
com <- com %>% mutate(cih90 = Effect + 1.645*se)
com <- com %>% mutate(cil90 = Effect - 1.645*se)

#bind results by treatment type
results_nid3 <- rbind(unlist(idt), unlist(cit), unlist(com))
colnames(results_nid3) <- c("Treatment_Type", "Effect", "se", "cih95", "cil95", "cih90", "cil90")
results_nid3 <- data.frame(results_nid3)
results_nid3$Effect <- as.numeric(results_nid3$Effect)
results_nid3$cil95 <- as.numeric(results_nid3$cil95)
results_nid3$cil90 <- as.numeric(results_nid3$cil90)
results_nid3$cih95 <- as.numeric(results_nid3$cih95)
results_nid3$cih90 <- as.numeric(results_nid3$cih90)
results_nid3$Quintile<- 3

#quintile 4

data1_nid4 <- subset(data1, data1$nat_quint==4)
data2_nid4 <- subset(data2, data2$nat_quint==4)
data3_nid4 <- subset(data3, data3$nat_quint==4)

#models by treatment type
lm1a <- lm_robust(natid_pca ~ as.factor(trt) + education + income + newjob, data=data1_nid4, se_type="HC3")
lm2a <- lm_robust(natid_pca ~ as.factor(trt) + lr_1 + natid_covar1_1 + income + newjob, data = data2_nid4, se_type="HC3")
lm3a <- lm_robust(natid_pca ~ as.factor(trt) + lr_1 + unemployed, data = data3_nid4, se_type="HC3")

idt <- c("Identity Threat\nvs. Control", lm1a$coefficients[2]) 
idt[2] <- lapply(idt[2], as.numeric); idt <- as.data.frame(idt)
colnames(idt) <- c("Treatment_Type", "Effect")
idt$se <- (lm1a$std.error)[2]
idt <- idt %>% mutate(cih95 = Effect + 1.96*se)
idt <- idt %>% mutate(cil95 = Effect - 1.96*se)
idt <- idt %>% mutate(cih90 = Effect + 1.645*se)
idt <- idt %>% mutate(cil90 = Effect - 1.645*se)

cit <- c("Whataboutism\nvs. Identity Threat", lm2a$coefficients[2]) 
cit[2] <- lapply(cit[2], as.numeric); cit <- as.data.frame(cit)
colnames(cit) <- c("Treatment_Type", "Effect")
cit$se <- (lm2a$std.error)[2]
cit <- cit %>% mutate(cih95 = Effect + 1.96*se)
cit <- cit %>% mutate(cil95 = Effect - 1.96*se)
cit <- cit %>% mutate(cih90 = Effect + 1.645*se)
cit <- cit %>% mutate(cil90 = Effect - 1.645*se)

com <- c("Whataboutism\nvs. Control", lm3a$coefficients[2]) 
com[2] <- lapply(com[2], as.numeric); com <- as.data.frame(com)
colnames(com) <- c("Treatment_Type", "Effect")
com$se <- (lm3a$std.error)[2]
com <- com %>% mutate(cih95 = Effect + 1.96*se)
com <- com %>% mutate(cil95 = Effect - 1.96*se)
com <- com %>% mutate(cih90 = Effect + 1.645*se)
com <- com %>% mutate(cil90 = Effect - 1.645*se)

#bind results by treatment type
results_nid4 <- rbind(unlist(idt), unlist(cit), unlist(com))
colnames(results_nid4) <- c("Treatment_Type", "Effect", "se", "cih95", "cil95", "cih90", "cil90")
results_nid4 <- data.frame(results_nid4)
results_nid4$Effect <- as.numeric(results_nid4$Effect)
results_nid4$cil95 <- as.numeric(results_nid4$cil95)
results_nid4$cil90 <- as.numeric(results_nid4$cil90)
results_nid4$cih95 <- as.numeric(results_nid4$cih95)
results_nid4$cih90 <- as.numeric(results_nid4$cih90)
results_nid4$Quintile<- 4

#quintile 5

data1_nid5 <- subset(data1, data1$nat_quint==5)
data2_nid5 <- subset(data2, data2$nat_quint==5)
data3_nid5 <- subset(data3, data3$nat_quint==5)

#models by treatment type
lm1a <- lm_robust(natid_pca ~ as.factor(trt) + education + income + newjob, data=data1_nid5, se_type="HC3")
lm2a <- lm_robust(natid_pca ~ as.factor(trt) + lr_1 + natid_covar1_1 + income + newjob, data = data2_nid5, se_type="HC3")
lm3a <- lm_robust(natid_pca ~ as.factor(trt) + lr_1 + unemployed, data = data3_nid5, se_type="HC3")

idt <- c("Identity Threat\nvs. Control", lm1a$coefficients[2]) 
idt[2] <- lapply(idt[2], as.numeric); idt <- as.data.frame(idt)
colnames(idt) <- c("Treatment_Type", "Effect")
idt$se <- (lm1a$std.error)[2]
idt <- idt %>% mutate(cih95 = Effect + 1.96*se)
idt <- idt %>% mutate(cil95 = Effect - 1.96*se)
idt <- idt %>% mutate(cih90 = Effect + 1.645*se)
idt <- idt %>% mutate(cil90 = Effect - 1.645*se)

cit <- c("Whataboutism\nvs. Identity Threat", lm2a$coefficients[2]) 
cit[2] <- lapply(cit[2], as.numeric); cit <- as.data.frame(cit)
colnames(cit) <- c("Treatment_Type", "Effect")
cit$se <- (lm2a$std.error)[2]
cit <- cit %>% mutate(cih95 = Effect + 1.96*se)
cit <- cit %>% mutate(cil95 = Effect - 1.96*se)
cit <- cit %>% mutate(cih90 = Effect + 1.645*se)
cit <- cit %>% mutate(cil90 = Effect - 1.645*se)

com <- c("Whataboutism\nvs. Control", lm3a$coefficients[2]) 
com[2] <- lapply(com[2], as.numeric); com <- as.data.frame(com)
colnames(com) <- c("Treatment_Type", "Effect")
com$se <- (lm3a$std.error)[2]
com <- com %>% mutate(cih95 = Effect + 1.96*se)
com <- com %>% mutate(cil95 = Effect - 1.96*se)
com <- com %>% mutate(cih90 = Effect + 1.645*se)
com <- com %>% mutate(cil90 = Effect - 1.645*se)

#bind results by treatment type
results_nid5 <- rbind(unlist(idt), unlist(cit), unlist(com))
colnames(results_nid5) <- c("Treatment_Type", "Effect", "se", "cih95", "cil95", "cih90", "cil90")
results_nid5 <- data.frame(results_nid5)
results_nid5$Effect <- as.numeric(results_nid5$Effect)
results_nid5$cil95 <- as.numeric(results_nid5$cil95)
results_nid5$cil90 <- as.numeric(results_nid5$cil90)
results_nid5$cih95 <- as.numeric(results_nid5$cih95)
results_nid5$cih90 <- as.numeric(results_nid5$cih90)
results_nid5$Quintile<- 5

#combine quintiles and plot
nid_hte <- rbind(results_nid1, results_nid2, results_nid3, results_nid4, results_nid5)
nid_hte$Effect <- as.numeric(as.character(nid_hte$Effect))

ni4 <- ggplot(nid_hte, aes(x = Quintile, y = Effect, colour=factor(Treatment_Type, level=factor_levels))) +
  geom_point(stat="identity", position=position_dodge(width=0.5), size=3) + 
  labs(tag = "D") +
  scale_y_continuous(limits = c(-0.55, 0.35)) +
  geom_linerange(aes(ymin = cil95, ymax = cih95), 
                 position = position_dodge(width=0.5,),
                 size = 1) +
  geom_linerange(aes(ymin = cil90, ymax = cih90), 
                 position = position_dodge(width=0.5),
                 size = 2)+
  ylab("Treatment Effect") +
  xlab("National Identification") +
  geom_hline(yintercept=0, linetype="dashed", color="blue") +
  labs(color='Update Type') +
  theme_bw()+
  theme(
    axis.text=element_text(size=16),
    plot.title = element_text(size=25),
    axis.title=element_text(size=18),
    legend.text=element_text(size=12),
    legend.title=element_text(size=16),
    legend.key.height = unit(1.5, 'cm')
  )

####Combined Panel for Figure D4####

(ni1 | ni2) /
  (ni3 | ni4)+
  plot_layout(guides = "collect") &
  theme(legend.position = "right")

ggsave("figD4.png", width = 12, height = 10, dpi = 300)

####TABLE D5 - HTEs for Perception of Empire####
####TABLE D5(A) - Perception of Empire by Left/Right Self-Placement####

#Left-wing
data1_left <- subset(data1, data1$lr_1 < 5)
data2_left <- subset(data2, data2$lr_1 < 5)
data3_left <- subset(data3, data3$lr_1 < 5)

#models by treatment type
lm1a <- lm_robust(empire_1 ~ as.factor(trt) + education + income + newjob, data=data1_left, se_type="HC3")
lm2a <- lm_robust(empire_1 ~ as.factor(trt) + lr_1 + natid_covar1_1 + income + newjob, data = data2_left, se_type="HC3")
lm3a <- lm_robust(empire_1 ~ as.factor(trt) + lr_1 + unemployed, data = data3_left, se_type="HC3")

idt <- c("Identity Threat\nvs. Control", lm1a$coefficients[2])
idt[2] <- lapply(idt[2], as.numeric); idt <- as.data.frame(idt)
colnames(idt) <- c("Treatment_Type", "Effect")
idt$se <- lm1a$std.error[2]
idt <- idt %>% 
  mutate(
    cih_95 = Effect + z_95 * se,
    cil_95 = Effect - z_95 * se,
    cih_90 = Effect + z_90 * se,
    cil_90 = Effect - z_90 * se
  )

cit <- c("Whataboutism\nvs. Identity Threat", lm2a$coefficients[2])
cit[2] <- lapply(cit[2], as.numeric); cit <- as.data.frame(cit)
colnames(cit) <- c("Treatment_Type", "Effect")
cit$se <- lm2a$std.error[2]
cit <- cit %>% 
  mutate(
    cih_95 = Effect + z_95 * se,
    cil_95 = Effect - z_95 * se,
    cih_90 = Effect + z_90 * se,
    cil_90 = Effect - z_90 * se
  )

com <- c("Whataboutism\nvs. Control", lm3a$coefficients[2])
com[2] <- lapply(com[2], as.numeric); com <- as.data.frame(com)
colnames(com) <- c("Treatment_Type", "Effect")
com$se <- lm3a$std.error[2]
com <- com %>% 
  mutate(
    cih_95 = Effect + z_95 * se,
    cil_95 = Effect - z_95 * se,
    cih_90 = Effect + z_90 * se,
    cil_90 = Effect - z_90 * se
  )

#bind results by treatment type
results_left <- rbind(unlist(idt), unlist(cit), unlist(com))
colnames(results_left) <- c("Treatment_Type", "Effect", "se", "cih_95", "cil_95", "cih_90", "cil_90")
results_left <- data.frame(results_left)
results_left$Effect <- as.numeric(results_left$Effect)
results_left$cil_95 <- as.numeric(results_left$cil_95)
results_left$cih_95 <- as.numeric(results_left$cih_95)
results_left$cil_90 <- as.numeric(results_left$cil_90)
results_left$cih_90 <- as.numeric(results_left$cih_90)
results_left$HTE <- "Left"

#Right-wing
data1_right <- subset(data1, data1$lr_1 >= 6)
data2_right <- subset(data2, data2$lr_1 >= 6)
data3_right <- subset(data3, data3$lr_1 >= 6)

#models by treatment type
lm1a <- lm_robust(empire_1 ~ as.factor(trt) + education + income + newjob, data=data1_right, se_type="HC3")
lm2a <- lm_robust(empire_1 ~ as.factor(trt) + lr_1 + natid_covar1_1 + income + newjob, data = data2_right, se_type="HC3")
lm3a <- lm_robust(empire_1 ~ as.factor(trt) + lr_1 + unemployed, data = data3_right, se_type="HC3")

idt <- c("Identity Threat\nvs. Control", lm1a$coefficients[2])
idt[2] <- lapply(idt[2], as.numeric); idt <- as.data.frame(idt)
colnames(idt) <- c("Treatment_Type", "Effect")
idt$se <- lm1a$std.error[2]
idt <- idt %>% 
  mutate(
    cih_95 = Effect + z_95 * se,
    cil_95 = Effect - z_95 * se,
    cih_90 = Effect + z_90 * se,
    cil_90 = Effect - z_90 * se
  )

cit <- c("Whataboutism\nvs. Identity Threat", lm2a$coefficients[2])
cit[2] <- lapply(cit[2], as.numeric); cit <- as.data.frame(cit)
colnames(cit) <- c("Treatment_Type", "Effect")
cit$se <- lm2a$std.error[2]
cit <- cit %>% 
  mutate(
    cih_95 = Effect + z_95 * se,
    cil_95 = Effect - z_95 * se,
    cih_90 = Effect + z_90 * se,
    cil_90 = Effect - z_90 * se
  )

com <- c("Whataboutism\nvs. Control", lm3a$coefficients[2])
com[2] <- lapply(com[2], as.numeric); com <- as.data.frame(com)
colnames(com) <- c("Treatment_Type", "Effect")
com$se <- lm3a$std.error[2]
com <- com %>% 
  mutate(
    cih_95 = Effect + z_95 * se,
    cil_95 = Effect - z_95 * se,
    cih_90 = Effect + z_90 * se,
    cil_90 = Effect - z_90 * se
  )

#bind results by treatment type
results_right <- rbind(unlist(idt), unlist(cit), unlist(com))
colnames(results_right) <- c("Treatment_Type", "Effect", "se", "cih_95", "cil_95", "cih_90", "cil_90")
results_right <- data.frame(results_right)
results_right$Effect <- as.numeric(results_right$Effect)
results_right$cil_95 <- as.numeric(results_right$cil_95)
results_right$cih_95 <- as.numeric(results_right$cih_95)
results_right$cil_90 <- as.numeric(results_right$cil_90)
results_right$cih_90 <- as.numeric(results_right$cih_90)
results_right$HTE <- "Right"

#bind left and right results and plot
lr_hte <- rbind(results_left, results_right)
lr_hte$HTE <- factor(lr_hte$HTE, levels = c("Left", "Right"))

em1 <- ggplot(lr_hte, aes(x = HTE, y = Effect, colour=factor(Treatment_Type, level=factor_levels))) +
  geom_point(stat="identity", position=position_dodge(width=0.3), size=3) + 
  labs(tag = "A") +
  scale_y_continuous(limits = c(-1.7, 1.3)) +
  geom_linerange(aes(ymin = cil_95, ymax = cih_95), 
                 position = position_dodge(width=0.3,),
                 size = 1) +
  geom_linerange(aes(ymin = cil_90, ymax = cih_90), 
                 position = position_dodge(width=0.3),
                 size = 2)+
  ylab("Treatment Effect") +
  xlab("Left-Right Self-Placement") +
  geom_hline(yintercept=0, linetype="dashed", color="blue") +
  labs(color='Update Type') +
  theme_bw()+
  theme(
    axis.text=element_text(size=16),
    plot.title = element_text(size=25),
    axis.title=element_text(size=18),
    legend.text=element_text(size=12),
    legend.title=element_text(size=16),
    legend.key.height = unit(1.5, 'cm')
  )

####TABLE D5(B) - Perception of Empire by Age####

#age group 1

data1_age1 <- subset(data1, data1$age<31)
data2_age1 <- subset(data2, data2$age<31)
data3_age1 <- subset(data3, data3$age<31)

#models by treatment type
lm1a <- lm_robust(empire_1 ~ as.factor(trt) + education + income + newjob, data=data1_age1, se_type="HC3")
lm2a <- lm_robust(empire_1 ~ as.factor(trt) + lr_1  + natid_covar1_1 + income + newjob, data=data2_age1, se_type="HC3")
lm3a <- lm_robust(empire_1 ~ as.factor(trt) + lr_1 + unemployed, data=data3_age1, se_type="HC3")

idt <- c("Identity Threat\nvs. Control", lm1a$coefficients[2]) 
idt[2] <- lapply(idt[2], as.numeric); idt <- as.data.frame(idt)
colnames(idt) <- c("Treatment_Type", "Effect")
idt$se <- (lm1a$std.error)[2]
idt <- idt %>% mutate(cih90 = Effect + 1.645*se)
idt <- idt %>% mutate(cil90 = Effect - 1.645*se)
idt <- idt %>% mutate(cih95 = Effect + 1.96*se)
idt <- idt %>% mutate(cil95 = Effect - 1.96*se)

cit <- c("Whataboutism\nvs. Identity Threat", lm2a$coefficients[2]) 
cit[2] <- lapply(cit[2], as.numeric); cit <- as.data.frame(cit)
colnames(cit) <- c("Treatment_Type", "Effect")
cit$se <- (lm2a$std.error)[2]
cit <- cit %>% mutate(cih90 = Effect + 1.645*se)
cit <- cit %>% mutate(cil90 = Effect - 1.645*se)
cit <- cit %>% mutate(cih95 = Effect + 1.96*se)
cit <- cit %>% mutate(cil95 = Effect - 1.96*se)

com <- c("Whataboutism\nvs. Control", lm3a$coefficients[2]) 
com[2] <- lapply(com[2], as.numeric); com <- as.data.frame(com)
colnames(com) <- c("Treatment_Type", "Effect")
com$se <- (lm3a$std.error)[2]
com <- com %>% mutate(cih90 = Effect + 1.645*se)
com <- com %>% mutate(cil90 = Effect - 1.645*se)
com <- com %>% mutate(cih95 = Effect + 1.96*se)
com <- com %>% mutate(cil95 = Effect - 1.96*se)

#bind results by treatment type
results_age1 <- rbind(unlist(idt), unlist(cit), unlist(com))
colnames(results_age1) <- c("Treatment_Type", "Effect", "se", "cih90", "cil90", "cih95", "cil95")
results_age1 <- data.frame(results_age1)
results_age1$Effect <- as.numeric(results_age1$Effect)
results_age1$cil90 <- as.numeric(results_age1$cil90)
results_age1$cih90 <- as.numeric(results_age1$cih90)
results_age1$cil95 <- as.numeric(results_age1$cil95)
results_age1$cih95 <- as.numeric(results_age1$cih95)

results_age1$Age<- "18-30"

#age group 2

data1_age2 <- subset(data1, data1$age>30 & data1$age<41)
data2_age2 <- subset(data2, data2$age>30 & data2$age<41)
data3_age2 <- subset(data3, data3$age>30 & data3$age<41)

#models by treatment type
lm1a <- lm_robust(empire_1 ~ as.factor(trt) + education + income + newjob, data=data1_age2, se_type="HC3")
lm2a <- lm_robust(empire_1 ~ as.factor(trt) + lr_1  + natid_covar1_1 + income + newjob, data=data2_age2, se_type="HC3")
lm3a <- lm_robust(empire_1 ~ as.factor(trt) + lr_1 + unemployed, data=data3_age2, se_type="HC3")

idt <- c("Identity Threat\nvs. Control", lm1a$coefficients[2]) 
idt[2] <- lapply(idt[2], as.numeric); idt <- as.data.frame(idt)
colnames(idt) <- c("Treatment_Type", "Effect")
idt$se <- (lm1a$std.error)[2]
idt <- idt %>% mutate(cih90 = Effect + 1.645*se)
idt <- idt %>% mutate(cil90 = Effect - 1.645*se)
idt <- idt %>% mutate(cih95 = Effect + 1.96*se)
idt <- idt %>% mutate(cil95 = Effect - 1.96*se)

cit <- c("Whataboutism\nvs. Identity Threat", lm2a$coefficients[2]) 
cit[2] <- lapply(cit[2], as.numeric); cit <- as.data.frame(cit)
colnames(cit) <- c("Treatment_Type", "Effect")
cit$se <- (lm2a$std.error)[2]
cit <- cit %>% mutate(cih90 = Effect + 1.645*se)
cit <- cit %>% mutate(cil90 = Effect - 1.645*se)
cit <- cit %>% mutate(cih95 = Effect + 1.96*se)
cit <- cit %>% mutate(cil95 = Effect - 1.96*se)

com <- c("Whataboutism\nvs. Control", lm3a$coefficients[2]) 
com[2] <- lapply(com[2], as.numeric); com <- as.data.frame(com)
colnames(com) <- c("Treatment_Type", "Effect")
com$se <- (lm3a$std.error)[2]
com <- com %>% mutate(cih90 = Effect + 1.645*se)
com <- com %>% mutate(cil90 = Effect - 1.645*se)
com <- com %>% mutate(cih95 = Effect + 1.96*se)
com <- com %>% mutate(cil95 = Effect - 1.96*se)

#bind results by treatment type
results_age2 <- rbind(unlist(idt), unlist(cit), unlist(com))
colnames(results_age2) <- c("Treatment_Type", "Effect", "se", "cih90", "cil90", "cih95", "cil95")
results_age2 <- data.frame(results_age2)
results_age2$Effect <- as.numeric(results_age2$Effect)
results_age2$cil90 <- as.numeric(results_age2$cil90)
results_age2$cih90 <- as.numeric(results_age2$cih90)
results_age2$cil95 <- as.numeric(results_age2$cil95)
results_age2$cih95 <- as.numeric(results_age2$cih95)
results_age2$Age<- "31-40"

#age group 3

data1_age3 <- subset(data1, data1$age>40 & data1$age<51)
data2_age3 <- subset(data2, data2$age>40 & data2$age<51)
data3_age3 <- subset(data3, data3$age>40 & data3$age<51)

#models by treatment type
lm1a <- lm_robust(empire_1 ~ as.factor(trt) + education + income + newjob, data=data1_age3, se_type="HC3")
lm2a <- lm_robust(empire_1 ~ as.factor(trt) + lr_1  + natid_covar1_1 + income + newjob, data=data2_age3, se_type="HC3")
lm3a <- lm_robust(empire_1 ~ as.factor(trt) + lr_1 + unemployed, data=data3_age3, se_type="HC3")

idt <- c("Identity Threat\nvs. Control", lm1a$coefficients[2]) 
idt[2] <- lapply(idt[2], as.numeric); idt <- as.data.frame(idt)
colnames(idt) <- c("Treatment_Type", "Effect")
idt$se <- (lm1a$std.error)[2]
idt <- idt %>% mutate(cih90 = Effect + 1.645*se)
idt <- idt %>% mutate(cil90 = Effect - 1.645*se)
idt <- idt %>% mutate(cih95 = Effect + 1.96*se)
idt <- idt %>% mutate(cil95 = Effect - 1.96*se)

cit <- c("Whataboutism\nvs. Identity Threat", lm2a$coefficients[2]) 
cit[2] <- lapply(cit[2], as.numeric); cit <- as.data.frame(cit)
colnames(cit) <- c("Treatment_Type", "Effect")
cit$se <- (lm2a$std.error)[2]
cit <- cit %>% mutate(cih90 = Effect + 1.645*se)
cit <- cit %>% mutate(cil90 = Effect - 1.645*se)
cit <- cit %>% mutate(cih95 = Effect + 1.96*se)
cit <- cit %>% mutate(cil95 = Effect - 1.96*se)

com <- c("Whataboutism\nvs. Control", lm3a$coefficients[2]) 
com[2] <- lapply(com[2], as.numeric); com <- as.data.frame(com)
colnames(com) <- c("Treatment_Type", "Effect")
com$se <- (lm3a$std.error)[2]
com <- com %>% mutate(cih90 = Effect + 1.645*se)
com <- com %>% mutate(cil90 = Effect - 1.645*se)
com <- com %>% mutate(cih95 = Effect + 1.96*se)
com <- com %>% mutate(cil95 = Effect - 1.96*se)

#bind results by treatment type
results_age3 <- rbind(unlist(idt), unlist(cit), unlist(com))
colnames(results_age3) <- c("Treatment_Type", "Effect", "se", "cih90", "cil90", "cih95", "cil95")
results_age3 <- data.frame(results_age3)
results_age3$Effect <- as.numeric(results_age3$Effect)
results_age3$cil90 <- as.numeric(results_age3$cil90)
results_age3$cih90 <- as.numeric(results_age3$cih90)
results_age3$cil95 <- as.numeric(results_age3$cil95)
results_age3$cih95 <- as.numeric(results_age3$cih95)
results_age3$Age<- "41-50"

#age group 4

data1_age4 <- subset(data1, data1$age>50 & data1$age<61)
data2_age4 <- subset(data2, data2$age>50 & data2$age<61)
data3_age4 <- subset(data3, data3$age>50 & data3$age<61)

#models by treatment type
lm1a <- lm_robust(empire_1 ~ as.factor(trt) + education + income + newjob, data=data1_age4, se_type="HC3")
lm2a <- lm_robust(empire_1 ~ as.factor(trt) + lr_1  + natid_covar1_1 + income + newjob, data=data2_age4, se_type="HC3")
lm3a <- lm_robust(empire_1 ~ as.factor(trt) + lr_1 + unemployed, data=data3_age4, se_type="HC3")

idt <- c("Identity Threat\nvs. Control", lm1a$coefficients[2]) 
idt[2] <- lapply(idt[2], as.numeric); idt <- as.data.frame(idt)
colnames(idt) <- c("Treatment_Type", "Effect")
idt$se <- (lm1a$std.error)[2]
idt <- idt %>% mutate(cih90 = Effect + 1.645*se)
idt <- idt %>% mutate(cil90 = Effect - 1.645*se)
idt <- idt %>% mutate(cih95 = Effect + 1.96*se)
idt <- idt %>% mutate(cil95 = Effect - 1.96*se)

cit <- c("Whataboutism\nvs. Identity Threat", lm2a$coefficients[2]) 
cit[2] <- lapply(cit[2], as.numeric); cit <- as.data.frame(cit)
colnames(cit) <- c("Treatment_Type", "Effect")
cit$se <- (lm2a$std.error)[2]
cit <- cit %>% mutate(cih90 = Effect + 1.645*se)
cit <- cit %>% mutate(cil90 = Effect - 1.645*se)
cit <- cit %>% mutate(cih95 = Effect + 1.96*se)
cit <- cit %>% mutate(cil95 = Effect - 1.96*se)

com <- c("Whataboutism\nvs. Control", lm3a$coefficients[2]) 
com[2] <- lapply(com[2], as.numeric); com <- as.data.frame(com)
colnames(com) <- c("Treatment_Type", "Effect")
com$se <- (lm3a$std.error)[2]
com <- com %>% mutate(cih90 = Effect + 1.645*se)
com <- com %>% mutate(cil90 = Effect - 1.645*se)
com <- com %>% mutate(cih95 = Effect + 1.96*se)
com <- com %>% mutate(cil95 = Effect - 1.96*se)

#bind results by treatment type
results_age4 <- rbind(unlist(idt), unlist(cit), unlist(com))
colnames(results_age4) <- c("Treatment_Type", "Effect", "se", "cih90", "cil90", "cih95", "cil95")
results_age4 <- data.frame(results_age4)
results_age4$Effect <- as.numeric(results_age4$Effect)
results_age4$cil90 <- as.numeric(results_age4$cil90)
results_age4$cih90 <- as.numeric(results_age4$cih90)
results_age4$cil95 <- as.numeric(results_age4$cil95)
results_age4$cih95 <- as.numeric(results_age4$cih95)
results_age4$Age<- "51-60"

#age group 5

data1_age5 <- subset(data1, data1$age>60)
data2_age5 <- subset(data2, data2$age>60)
data3_age5 <- subset(data3, data3$age>60)

#models by treatment type
lm1a <- lm_robust(empire_1 ~ as.factor(trt) + education + income + newjob, data=data1_age5, se_type="HC3")
lm2a <- lm_robust(empire_1 ~ as.factor(trt) + lr_1  + natid_covar1_1 + income + newjob, data=data2_age5, se_type="HC3")
lm3a <- lm_robust(empire_1 ~ as.factor(trt) + lr_1 + unemployed, data=data3_age5, se_type="HC3")

idt <- c("Identity Threat\nvs. Control", lm1a$coefficients[2]) 
idt[2] <- lapply(idt[2], as.numeric); idt <- as.data.frame(idt)
colnames(idt) <- c("Treatment_Type", "Effect")
idt$se <- (lm1a$std.error)[2]
idt <- idt %>% mutate(cih90 = Effect + 1.645*se)
idt <- idt %>% mutate(cil90 = Effect - 1.645*se)
idt <- idt %>% mutate(cih95 = Effect + 1.96*se)
idt <- idt %>% mutate(cil95 = Effect - 1.96*se)

cit <- c("Whataboutism\nvs. Identity Threat", lm2a$coefficients[2]) 
cit[2] <- lapply(cit[2], as.numeric); cit <- as.data.frame(cit)
colnames(cit) <- c("Treatment_Type", "Effect")
cit$se <- (lm2a$std.error)[2]
cit <- cit %>% mutate(cih90 = Effect + 1.645*se)
cit <- cit %>% mutate(cil90 = Effect - 1.645*se)
cit <- cit %>% mutate(cih95 = Effect + 1.96*se)
cit <- cit %>% mutate(cil95 = Effect - 1.96*se)

com <- c("Whataboutism\nvs. Control", lm3a$coefficients[2]) 
com[2] <- lapply(com[2], as.numeric); com <- as.data.frame(com)
colnames(com) <- c("Treatment_Type", "Effect")
com$se <- (lm3a$std.error)[2]
com <- com %>% mutate(cih90 = Effect + 1.645*se)
com <- com %>% mutate(cil90 = Effect - 1.645*se)
com <- com %>% mutate(cih95 = Effect + 1.96*se)
com <- com %>% mutate(cil95 = Effect - 1.96*se)

#bind results by treatment type
results_age5 <- rbind(unlist(idt), unlist(cit), unlist(com))
colnames(results_age5) <- c("Treatment_Type", "Effect", "se", "cih90", "cil90", "cih95", "cil95")
results_age5 <- data.frame(results_age5)
results_age5$Effect <- as.numeric(results_age5$Effect)
results_age5$cil90 <- as.numeric(results_age5$cil90)
results_age5$cih90 <- as.numeric(results_age5$cih90)
results_age5$cil95 <- as.numeric(results_age5$cil95)
results_age5$cih95 <- as.numeric(results_age5$cih95)
results_age5$Age<- "61+"

#combine age groups and plot

age_hte <- rbind(results_age1, results_age2, results_age3, results_age4, results_age5)
age_hte$Age <- factor(age_hte$Age, levels = c("18-30", "31-40", "41-50", "51-60", "61+"))

em2 <- ggplot(age_hte, aes(x = Age, y = Effect, colour=factor(Treatment_Type, level=factor_levels))) +
  geom_point(stat="identity", position=position_dodge(width=0.5), size=3) + 
  labs(tag = "B") +
  scale_y_continuous(limits = c(-1.7, 1.3)) +
  geom_linerange(aes(ymin = cil95, ymax = cih95), 
                 position = position_dodge(width=0.5,),
                 size = 1) +
  geom_linerange(aes(ymin = cil90, ymax = cih90), 
                 position = position_dodge(width=0.5),
                 size = 2) +
  ylab("Treatment Effect") +
  xlab("Age") +
  geom_hline(yintercept=0, linetype="dashed", color="blue") +
  labs(color='Update Type') +
  theme_bw()+
  theme(
    axis.text=element_text(size=16),
    plot.title = element_text(size=25),
    axis.title=element_text(size=18),
    legend.text=element_text(size=12),
    legend.title=element_text(size=16),
    legend.key.height = unit(1.5, 'cm')
  )

####TABLE D5(C) - Perception of Empire by Gender####

# women
data1_f <- subset(data1, data1$female == 1)
data2_f <- subset(data2, data2$female == 1)
data3_f <- subset(data3, data3$female == 1)

#models by treatment type
lm1a <- lm_robust(empire_1 ~ as.factor(trt) + education + income + newjob, data=data1_f, se_type="HC3")
lm2a <- lm_robust(empire_1 ~ as.factor(trt) + lr_1 + natid_covar1_1 + income + newjob, data = data2_f, se_type="HC3")
lm3a <- lm_robust(empire_1 ~ as.factor(trt) + lr_1 + unemployed, data = data3_f, se_type="HC3")

idt <- data.frame(
  Treatment_Type = "Identity Threat\nvs. Control",
  Effect = lm1a$coefficients[2],
  se = lm1a$std.error[2]
)
idt <- idt %>% 
  mutate(
    ci90h = Effect + 1.645 * se,
    ci90l = Effect - 1.645 * se,
    ci95h = Effect + 1.96 * se,
    ci95l = Effect - 1.96 * se
  )

cit <- data.frame(
  Treatment_Type = "Whataboutism\nvs. Identity Threat",
  Effect = lm2a$coefficients[2],
  se = lm2a$std.error[2]
)
cit <- cit %>% 
  mutate(
    ci90h = Effect + 1.645 * se,
    ci90l = Effect - 1.645 * se,
    ci95h = Effect + 1.96 * se,
    ci95l = Effect - 1.96 * se
  )

com <- data.frame(
  Treatment_Type = "Whataboutism\nvs. Control",
  Effect = lm3a$coefficients[2],
  se = lm3a$std.error[2]
)
com <- com %>% 
  mutate(
    ci90h = Effect + 1.645 * se,
    ci90l = Effect - 1.645 * se,
    ci95h = Effect + 1.96 * se,
    ci95l = Effect - 1.96 * se
  )

#bind results by treatment type
results_f <- rbind(idt, cit, com)
results_f$Gender <- "Female"

# men
data1_m <- subset(data1, data1$female == 0)
data2_m <- subset(data2, data2$female == 0)
data3_m <- subset(data3, data3$female == 0)

#models by treatment type
lm1a <- lm_robust(empire_1 ~ as.factor(trt) + education + income + newjob, data=data1_m, se_type="HC3")
lm2a <- lm_robust(empire_1 ~ as.factor(trt) + lr_1 + natid_covar1_1 + income + newjob, data = data2_m, se_type="HC3")
lm3a <- lm_robust(empire_1 ~ as.factor(trt) + lr_1 + unemployed, data = data3_m, se_type="HC3")

idt <- data.frame(
  Treatment_Type = "Identity Threat\nvs. Control",
  Effect = lm1a$coefficients[2],
  se = lm1a$std.error[2]
)
idt <- idt %>% 
  mutate(
    ci90h = Effect + 1.645 * se,
    ci90l = Effect - 1.645 * se,
    ci95h = Effect + 1.96 * se,
    ci95l = Effect - 1.96 * se
  )

cit <- data.frame(
  Treatment_Type = "Whataboutism\nvs. Identity Threat",
  Effect = lm2a$coefficients[2],
  se = lm2a$std.error[2]
)
cit <- cit %>% 
  mutate(
    ci90h = Effect + 1.645 * se,
    ci90l = Effect - 1.645 * se,
    ci95h = Effect + 1.96 * se,
    ci95l = Effect - 1.96 * se
  )

com <- data.frame(
  Treatment_Type = "Whataboutism\nvs. Control",
  Effect = lm3a$coefficients[2],
  se = lm3a$std.error[2]
)
com <- com %>% 
  mutate(
    ci90h = Effect + 1.645 * se,
    ci90l = Effect - 1.645 * se,
    ci95h = Effect + 1.96 * se,
    ci95l = Effect - 1.96 * se
  )

#bind results by treatment type
results_m <- rbind(idt, cit, com)
results_m$Gender <- "Male"

# Combine women/men and plot
gender_hte <- rbind(results_f, results_m)
gender_hte$Gender <- factor(gender_hte$Gender, levels = c("Female", "Male"))

em3 <- ggplot(gender_hte, aes(x = Gender, y = Effect, colour=factor(Treatment_Type, level=factor_levels))) +
  geom_point(stat="identity", position=position_dodge(width=0.3), size=3) + 
  labs(tag = "C") +
  scale_y_continuous(limits = c(-1.7, 1.3)) +
  geom_linerange(aes(ymin = ci95l, ymax = ci95h), 
                 position = position_dodge(width=0.3,),
                 size = 1) +
  geom_linerange(aes(ymin = ci90l, ymax = ci90h), 
                 position = position_dodge(width=0.3),
                 size = 2)+
  ylab("Treatment Effect") +
  xlab("Gender") +
  geom_hline(yintercept=0, linetype="dashed", color="blue") +
  labs(color='Update Type') +
  theme_bw()+
  theme(
    axis.text=element_text(size=16),
    plot.title = element_text(size=25),
    axis.title=element_text(size=18),
    legend.text=element_text(size=12),
    legend.title=element_text(size=16),
    legend.key.height = unit(1.5, 'cm')
  )


####TABLE D5(D) - Perception of Empire by Pre-Treatment National Identification####

#quintile 1

data1_nid1 <- subset(data1, data1$nat_quint==1)
data2_nid1 <- subset(data2, data2$nat_quint==1)
data3_nid1 <- subset(data3, data3$nat_quint==1)

#models by treatment type
lm1a <- lm_robust(empire_1 ~ as.factor(trt) + education + income + newjob, data=data1_nid1, se_type="HC3")
lm2a <- lm_robust(empire_1 ~ as.factor(trt) + lr_1 + natid_covar1_1 + income + newjob, data = data2_nid1, se_type="HC3")
lm3a <- lm_robust(empire_1 ~ as.factor(trt) + lr_1 + unemployed, data = data3_nid1, se_type="HC3")

idt <- c("Identity Threat\nvs. Control", lm1a$coefficients[2]) 
idt[2] <- lapply(idt[2], as.numeric); idt <- as.data.frame(idt)
colnames(idt) <- c("Treatment_Type", "Effect")
idt$se <- (lm1a$std.error)[2]
idt <- idt %>% mutate(cih95 = Effect + 1.96*se)
idt <- idt %>% mutate(cil95 = Effect - 1.96*se)
idt <- idt %>% mutate(cih90 = Effect + 1.645*se)
idt <- idt %>% mutate(cil90 = Effect - 1.645*se)

cit <- c("Whataboutism\nvs. Identity Threat", lm2a$coefficients[2]) 
cit[2] <- lapply(cit[2], as.numeric); cit <- as.data.frame(cit)
colnames(cit) <- c("Treatment_Type", "Effect")
cit$se <- (lm2a$std.error)[2]
cit <- cit %>% mutate(cih95 = Effect + 1.96*se)
cit <- cit %>% mutate(cil95 = Effect - 1.96*se)
cit <- cit %>% mutate(cih90 = Effect + 1.645*se)
cit <- cit %>% mutate(cil90 = Effect - 1.645*se)

com <- c("Whataboutism\nvs. Control", lm3a$coefficients[2]) 
com[2] <- lapply(com[2], as.numeric); com <- as.data.frame(com)
colnames(com) <- c("Treatment_Type", "Effect")
com$se <- (lm3a$std.error)[2]
com <- com %>% mutate(cih95 = Effect + 1.96*se)
com <- com %>% mutate(cil95 = Effect - 1.96*se)
com <- com %>% mutate(cih90 = Effect + 1.645*se)
com <- com %>% mutate(cil90 = Effect - 1.645*se)

#bind results by treatment type
results_nid1 <- rbind(unlist(idt), unlist(cit), unlist(com))
colnames(results_nid1) <- c("Treatment_Type", "Effect", "se", "cih95", "cil95", "cih90", "cil90")
results_nid1 <- data.frame(results_nid1)
results_nid1$Effect <- as.numeric(results_nid1$Effect)
results_nid1$cil95 <- as.numeric(results_nid1$cil95)
results_nid1$cil90 <- as.numeric(results_nid1$cil90)
results_nid1$cih95 <- as.numeric(results_nid1$cih95)
results_nid1$cih90 <- as.numeric(results_nid1$cih90)
results_nid1$Quintile<- 1

#quintile 2

data1_nid2 <- subset(data1, data1$nat_quint==2)
data2_nid2 <- subset(data2, data2$nat_quint==2)
data3_nid2 <- subset(data3, data3$nat_quint==2)

#models by treatment type
lm1a <- lm_robust(empire_1 ~ as.factor(trt) + education + income + newjob, data=data1_nid2, se_type="HC3")
lm2a <- lm_robust(empire_1 ~ as.factor(trt) + lr_1 + natid_covar1_1 + income + newjob, data = data2_nid2, se_type="HC3")
lm3a <- lm_robust(empire_1 ~ as.factor(trt) + lr_1 + unemployed, data = data3_nid2, se_type="HC3")

idt <- c("Identity Threat\nvs. Control", lm1a$coefficients[2]) 
idt[2] <- lapply(idt[2], as.numeric); idt <- as.data.frame(idt)
colnames(idt) <- c("Treatment_Type", "Effect")
idt$se <- (lm1a$std.error)[2]
idt <- idt %>% mutate(cih95 = Effect + 1.96*se)
idt <- idt %>% mutate(cil95 = Effect - 1.96*se)
idt <- idt %>% mutate(cih90 = Effect + 1.645*se)
idt <- idt %>% mutate(cil90 = Effect - 1.645*se)

cit <- c("Whataboutism\nvs. Identity Threat", lm2a$coefficients[2]) 
cit[2] <- lapply(cit[2], as.numeric); cit <- as.data.frame(cit)
colnames(cit) <- c("Treatment_Type", "Effect")
cit$se <- (lm2a$std.error)[2]
cit <- cit %>% mutate(cih95 = Effect + 1.96*se)
cit <- cit %>% mutate(cil95 = Effect - 1.96*se)
cit <- cit %>% mutate(cih90 = Effect + 1.645*se)
cit <- cit %>% mutate(cil90 = Effect - 1.645*se)

com <- c("Whataboutism\nvs. Control", lm3a$coefficients[2]) 
com[2] <- lapply(com[2], as.numeric); com <- as.data.frame(com)
colnames(com) <- c("Treatment_Type", "Effect")
com$se <- (lm3a$std.error)[2]
com <- com %>% mutate(cih95 = Effect + 1.96*se)
com <- com %>% mutate(cil95 = Effect - 1.96*se)
com <- com %>% mutate(cih90 = Effect + 1.645*se)
com <- com %>% mutate(cil90 = Effect - 1.645*se)

#bind results by treatment type
results_nid2 <- rbind(unlist(idt), unlist(cit), unlist(com))
colnames(results_nid2) <- c("Treatment_Type", "Effect", "se", "cih95", "cil95", "cih90", "cil90")
results_nid2 <- data.frame(results_nid2)
results_nid2$Effect <- as.numeric(results_nid2$Effect)
results_nid2$cil95 <- as.numeric(results_nid2$cil95)
results_nid2$cil90 <- as.numeric(results_nid2$cil90)
results_nid2$cih95 <- as.numeric(results_nid2$cih95)
results_nid2$cih90 <- as.numeric(results_nid2$cih90)
results_nid2$Quintile<- 2

#quintile 3

data1_nid3 <- subset(data1, data1$nat_quint==3)
data2_nid3 <- subset(data2, data2$nat_quint==3)
data3_nid3 <- subset(data3, data3$nat_quint==3)

#models by treatment type
lm1a <- lm_robust(empire_1 ~ as.factor(trt) + education + income + newjob, data=data1_nid3, se_type="HC3")
lm2a <- lm_robust(empire_1 ~ as.factor(trt) + lr_1 + natid_covar1_1 + income + newjob, data = data2_nid3, se_type="HC3")
lm3a <- lm_robust(empire_1 ~ as.factor(trt) + lr_1 + unemployed, data = data3_nid3, se_type="HC3")

idt <- c("Identity Threat\nvs. Control", lm1a$coefficients[2]) 
idt[2] <- lapply(idt[2], as.numeric); idt <- as.data.frame(idt)
colnames(idt) <- c("Treatment_Type", "Effect")
idt$se <- (lm1a$std.error)[2]
idt <- idt %>% mutate(cih95 = Effect + 1.96*se)
idt <- idt %>% mutate(cil95 = Effect - 1.96*se)
idt <- idt %>% mutate(cih90 = Effect + 1.645*se)
idt <- idt %>% mutate(cil90 = Effect - 1.645*se)

cit <- c("Whataboutism\nvs. Identity Threat", lm2a$coefficients[2]) 
cit[2] <- lapply(cit[2], as.numeric); cit <- as.data.frame(cit)
colnames(cit) <- c("Treatment_Type", "Effect")
cit$se <- (lm2a$std.error)[2]
cit <- cit %>% mutate(cih95 = Effect + 1.96*se)
cit <- cit %>% mutate(cil95 = Effect - 1.96*se)
cit <- cit %>% mutate(cih90 = Effect + 1.645*se)
cit <- cit %>% mutate(cil90 = Effect - 1.645*se)

com <- c("Whataboutism\nvs. Control", lm3a$coefficients[2]) 
com[2] <- lapply(com[2], as.numeric); com <- as.data.frame(com)
colnames(com) <- c("Treatment_Type", "Effect")
com$se <- (lm3a$std.error)[2]
com <- com %>% mutate(cih95 = Effect + 1.96*se)
com <- com %>% mutate(cil95 = Effect - 1.96*se)
com <- com %>% mutate(cih90 = Effect + 1.645*se)
com <- com %>% mutate(cil90 = Effect - 1.645*se)

#bind results by treatment type
results_nid3 <- rbind(unlist(idt), unlist(cit), unlist(com))
colnames(results_nid3) <- c("Treatment_Type", "Effect", "se", "cih95", "cil95", "cih90", "cil90")
results_nid3 <- data.frame(results_nid3)
results_nid3$Effect <- as.numeric(results_nid3$Effect)
results_nid3$cil95 <- as.numeric(results_nid3$cil95)
results_nid3$cil90 <- as.numeric(results_nid3$cil90)
results_nid3$cih95 <- as.numeric(results_nid3$cih95)
results_nid3$cih90 <- as.numeric(results_nid3$cih90)
results_nid3$Quintile<- 3

#quintile 4

data1_nid4 <- subset(data1, data1$nat_quint==4)
data2_nid4 <- subset(data2, data2$nat_quint==4)
data3_nid4 <- subset(data3, data3$nat_quint==4)

#models by treatment type
lm1a <- lm_robust(empire_1 ~ as.factor(trt) + education + income + newjob, data=data1_nid4, se_type="HC3")
lm2a <- lm_robust(empire_1 ~ as.factor(trt) + lr_1 + natid_covar1_1 + income + newjob, data = data2_nid4, se_type="HC3")
lm3a <- lm_robust(empire_1 ~ as.factor(trt) + lr_1 + unemployed, data = data3_nid4, se_type="HC3")

idt <- c("Identity Threat\nvs. Control", lm1a$coefficients[2]) 
idt[2] <- lapply(idt[2], as.numeric); idt <- as.data.frame(idt)
colnames(idt) <- c("Treatment_Type", "Effect")
idt$se <- (lm1a$std.error)[2]
idt <- idt %>% mutate(cih95 = Effect + 1.96*se)
idt <- idt %>% mutate(cil95 = Effect - 1.96*se)
idt <- idt %>% mutate(cih90 = Effect + 1.645*se)
idt <- idt %>% mutate(cil90 = Effect - 1.645*se)

cit <- c("Whataboutism\nvs. Identity Threat", lm2a$coefficients[2]) 
cit[2] <- lapply(cit[2], as.numeric); cit <- as.data.frame(cit)
colnames(cit) <- c("Treatment_Type", "Effect")
cit$se <- (lm2a$std.error)[2]
cit <- cit %>% mutate(cih95 = Effect + 1.96*se)
cit <- cit %>% mutate(cil95 = Effect - 1.96*se)
cit <- cit %>% mutate(cih90 = Effect + 1.645*se)
cit <- cit %>% mutate(cil90 = Effect - 1.645*se)

com <- c("Whataboutism\nvs. Control", lm3a$coefficients[2]) 
com[2] <- lapply(com[2], as.numeric); com <- as.data.frame(com)
colnames(com) <- c("Treatment_Type", "Effect")
com$se <- (lm3a$std.error)[2]
com <- com %>% mutate(cih95 = Effect + 1.96*se)
com <- com %>% mutate(cil95 = Effect - 1.96*se)
com <- com %>% mutate(cih90 = Effect + 1.645*se)
com <- com %>% mutate(cil90 = Effect - 1.645*se)

#bind results by treatment type
results_nid4 <- rbind(unlist(idt), unlist(cit), unlist(com))
colnames(results_nid4) <- c("Treatment_Type", "Effect", "se", "cih95", "cil95", "cih90", "cil90")
results_nid4 <- data.frame(results_nid4)
results_nid4$Effect <- as.numeric(results_nid4$Effect)
results_nid4$cil95 <- as.numeric(results_nid4$cil95)
results_nid4$cil90 <- as.numeric(results_nid4$cil90)
results_nid4$cih95 <- as.numeric(results_nid4$cih95)
results_nid4$cih90 <- as.numeric(results_nid4$cih90)
results_nid4$Quintile<- 4

#quintile 5

data1_nid5 <- subset(data1, data1$nat_quint==5)
data2_nid5 <- subset(data2, data2$nat_quint==5)
data3_nid5 <- subset(data3, data3$nat_quint==5)

#models by treatment type
lm1a <- lm_robust(empire_1 ~ as.factor(trt) + education + income + newjob, data=data1_nid5, se_type="HC3")
lm2a <- lm_robust(empire_1 ~ as.factor(trt) + lr_1 + natid_covar1_1 + income + newjob, data = data2_nid5, se_type="HC3")
lm3a <- lm_robust(empire_1 ~ as.factor(trt) + lr_1 + unemployed, data = data3_nid5, se_type="HC3")

idt <- c("Identity Threat\nvs. Control", lm1a$coefficients[2]) 
idt[2] <- lapply(idt[2], as.numeric); idt <- as.data.frame(idt)
colnames(idt) <- c("Treatment_Type", "Effect")
idt$se <- (lm1a$std.error)[2]
idt <- idt %>% mutate(cih95 = Effect + 1.96*se)
idt <- idt %>% mutate(cil95 = Effect - 1.96*se)
idt <- idt %>% mutate(cih90 = Effect + 1.645*se)
idt <- idt %>% mutate(cil90 = Effect - 1.645*se)

cit <- c("Whataboutism\nvs. Identity Threat", lm2a$coefficients[2]) 
cit[2] <- lapply(cit[2], as.numeric); cit <- as.data.frame(cit)
colnames(cit) <- c("Treatment_Type", "Effect")
cit$se <- (lm2a$std.error)[2]
cit <- cit %>% mutate(cih95 = Effect + 1.96*se)
cit <- cit %>% mutate(cil95 = Effect - 1.96*se)
cit <- cit %>% mutate(cih90 = Effect + 1.645*se)
cit <- cit %>% mutate(cil90 = Effect - 1.645*se)

com <- c("Whataboutism\nvs. Control", lm3a$coefficients[2]) 
com[2] <- lapply(com[2], as.numeric); com <- as.data.frame(com)
colnames(com) <- c("Treatment_Type", "Effect")
com$se <- (lm3a$std.error)[2]
com <- com %>% mutate(cih95 = Effect + 1.96*se)
com <- com %>% mutate(cil95 = Effect - 1.96*se)
com <- com %>% mutate(cih90 = Effect + 1.645*se)
com <- com %>% mutate(cil90 = Effect - 1.645*se)

#bind results by treatment type
results_nid5 <- rbind(unlist(idt), unlist(cit), unlist(com))
colnames(results_nid5) <- c("Treatment_Type", "Effect", "se", "cih95", "cil95", "cih90", "cil90")
results_nid5 <- data.frame(results_nid5)
results_nid5$Effect <- as.numeric(results_nid5$Effect)
results_nid5$cil95 <- as.numeric(results_nid5$cil95)
results_nid5$cil90 <- as.numeric(results_nid5$cil90)
results_nid5$cih95 <- as.numeric(results_nid5$cih95)
results_nid5$cih90 <- as.numeric(results_nid5$cih90)
results_nid5$Quintile<- 5

#combine quintiles and plot
nid_hte <- rbind(results_nid1, results_nid2, results_nid3, results_nid4, results_nid5)
nid_hte$Effect <- as.numeric(as.character(nid_hte$Effect))

em4 <- ggplot(nid_hte, aes(x = Quintile, y = Effect, colour=factor(Treatment_Type, level=factor_levels))) +
  geom_point(stat="identity", position=position_dodge(width=0.5), size=3) + 
  labs(tag = "D") +
  scale_y_continuous(limits = c(-1.7, 1.3)) +
  geom_linerange(aes(ymin = cil95, ymax = cih95), 
                 position = position_dodge(width=0.5,),
                 size = 1) +
  geom_linerange(aes(ymin = cil90, ymax = cih90), 
                 position = position_dodge(width=0.5),
                 size = 2)+
  ylab("Treatment Effect") +
  xlab("National Identification") +
  geom_hline(yintercept=0, linetype="dashed", color="blue") +
  labs(color='Update Type') +
  theme_bw()+
  theme(
    axis.text=element_text(size=16),
    plot.title = element_text(size=25),
    axis.title=element_text(size=18),
    legend.text=element_text(size=12),
    legend.title=element_text(size=16),
    legend.key.height = unit(1.5, 'cm')
  )

####Combined Panel for Figure D5####

(em1 | em2) /
  (em3 | em4)+
  plot_layout(guides = "collect") &
  theme(legend.position = "right")

ggsave("figD5.png", width = 12, height = 10, dpi = 300)

####TABLE D6 - HTEs for Whataboutism####
####TABLE D6(A) - Whataboutism by Left/Right Self-Placement####

#Left-wing
data1_left <- subset(data1, data1$lr_1 < 5)
data2_left <- subset(data2, data2$lr_1 < 5)
data3_left <- subset(data3, data3$lr_1 < 5)

#models by treatment type
lm1a <- lm_robust(comparison_1 ~ as.factor(trt) + education + income + newjob, data=data1_left, se_type="HC3")
lm2a <- lm_robust(comparison_1 ~ as.factor(trt) + lr_1 + natid_covar1_1 + income + newjob, data = data2_left, se_type="HC3")
lm3a <- lm_robust(comparison_1 ~ as.factor(trt) + lr_1 + unemployed, data = data3_left, se_type="HC3")

idt <- c("Identity Threat\nvs. Control", lm1a$coefficients[2])
idt[2] <- lapply(idt[2], as.numeric); idt <- as.data.frame(idt)
colnames(idt) <- c("Treatment_Type", "Effect")
idt$se <- lm1a$std.error[2]
idt <- idt %>% 
  mutate(
    cih_95 = Effect + z_95 * se,
    cil_95 = Effect - z_95 * se,
    cih_90 = Effect + z_90 * se,
    cil_90 = Effect - z_90 * se
  )

cit <- c("Whataboutism\nvs. Identity Threat", lm2a$coefficients[2])
cit[2] <- lapply(cit[2], as.numeric); cit <- as.data.frame(cit)
colnames(cit) <- c("Treatment_Type", "Effect")
cit$se <- lm2a$std.error[2]
cit <- cit %>% 
  mutate(
    cih_95 = Effect + z_95 * se,
    cil_95 = Effect - z_95 * se,
    cih_90 = Effect + z_90 * se,
    cil_90 = Effect - z_90 * se
  )

com <- c("Whataboutism\nvs. Control", lm3a$coefficients[2])
com[2] <- lapply(com[2], as.numeric); com <- as.data.frame(com)
colnames(com) <- c("Treatment_Type", "Effect")
com$se <- lm3a$std.error[2]
com <- com %>% 
  mutate(
    cih_95 = Effect + z_95 * se,
    cil_95 = Effect - z_95 * se,
    cih_90 = Effect + z_90 * se,
    cil_90 = Effect - z_90 * se
  )

#bind results by treatment type
results_left <- rbind(unlist(idt), unlist(cit), unlist(com))
colnames(results_left) <- c("Treatment_Type", "Effect", "se", "cih_95", "cil_95", "cih_90", "cil_90")
results_left <- data.frame(results_left)
results_left$Effect <- as.numeric(results_left$Effect)
results_left$cil_95 <- as.numeric(results_left$cil_95)
results_left$cih_95 <- as.numeric(results_left$cih_95)
results_left$cil_90 <- as.numeric(results_left$cil_90)
results_left$cih_90 <- as.numeric(results_left$cih_90)
results_left$HTE <- "Left"

#Right-wing
data1_right <- subset(data1, data1$lr_1 >= 6)
data2_right <- subset(data2, data2$lr_1 >= 6)
data3_right <- subset(data3, data3$lr_1 >= 6)

#models by treatment type
lm1a <- lm_robust(comparison_1 ~ as.factor(trt) + education + income + newjob, data=data1_right, se_type="HC3")
lm2a <- lm_robust(comparison_1 ~ as.factor(trt) + lr_1 + natid_covar1_1 + income + newjob, data = data2_right, se_type="HC3")
lm3a <- lm_robust(comparison_1 ~ as.factor(trt) + lr_1 + unemployed, data = data3_right, se_type="HC3")

idt <- c("Identity Threat\nvs. Control", lm1a$coefficients[2])
idt[2] <- lapply(idt[2], as.numeric); idt <- as.data.frame(idt)
colnames(idt) <- c("Treatment_Type", "Effect")
idt$se <- lm1a$std.error[2]
idt <- idt %>% 
  mutate(
    cih_95 = Effect + z_95 * se,
    cil_95 = Effect - z_95 * se,
    cih_90 = Effect + z_90 * se,
    cil_90 = Effect - z_90 * se
  )

cit <- c("Whataboutism\nvs. Identity Threat", lm2a$coefficients[2])
cit[2] <- lapply(cit[2], as.numeric); cit <- as.data.frame(cit)
colnames(cit) <- c("Treatment_Type", "Effect")
cit$se <- lm2a$std.error[2]
cit <- cit %>% 
  mutate(
    cih_95 = Effect + z_95 * se,
    cil_95 = Effect - z_95 * se,
    cih_90 = Effect + z_90 * se,
    cil_90 = Effect - z_90 * se
  )

com <- c("Whataboutism\nvs. Control", lm3a$coefficients[2])
com[2] <- lapply(com[2], as.numeric); com <- as.data.frame(com)
colnames(com) <- c("Treatment_Type", "Effect")
com$se <- lm3a$std.error[2]
com <- com %>% 
  mutate(
    cih_95 = Effect + z_95 * se,
    cil_95 = Effect - z_95 * se,
    cih_90 = Effect + z_90 * se,
    cil_90 = Effect - z_90 * se
  )

#bind results by treatment type
results_right <- rbind(unlist(idt), unlist(cit), unlist(com))
colnames(results_right) <- c("Treatment_Type", "Effect", "se", "cih_95", "cil_95", "cih_90", "cil_90")
results_right <- data.frame(results_right)
results_right$Effect <- as.numeric(results_right$Effect)
results_right$cil_95 <- as.numeric(results_right$cil_95)
results_right$cih_95 <- as.numeric(results_right$cih_95)
results_right$cil_90 <- as.numeric(results_right$cil_90)
results_right$cih_90 <- as.numeric(results_right$cih_90)
results_right$HTE <- "Right"

#bind left and right results and plot
lr_hte <- rbind(results_left, results_right)
lr_hte$HTE <- factor(lr_hte$HTE, levels = c("Left", "Right"))

wb1 <- ggplot(lr_hte, aes(x = HTE, y = Effect, colour=factor(Treatment_Type, level=factor_levels))) +
  geom_point(stat="identity", position=position_dodge(width=0.3), size=3) + 
  labs(tag = "A") +
  scale_y_continuous(limits = c(-1.55, 1.55)) +
  geom_linerange(aes(ymin = cil_95, ymax = cih_95), 
                 position = position_dodge(width=0.3,),
                 size = 1) +
  geom_linerange(aes(ymin = cil_90, ymax = cih_90), 
                 position = position_dodge(width=0.3),
                 size = 2)+
  ylab("Treatment Effect") +
  xlab("Left-Right Self-Placement") +
  geom_hline(yintercept=0, linetype="dashed", color="blue") +
  labs(color='Update Type') +
  theme_bw()+
  theme(
    axis.text=element_text(size=16),
    plot.title = element_text(size=25),
    axis.title=element_text(size=18),
    legend.text=element_text(size=12),
    legend.title=element_text(size=16),
    legend.key.height = unit(1.5, 'cm')
  )

####TABLE D6(B) - Whataboutism by Age####

#age group 1

data1_age1 <- subset(data1, data1$age<31)
data2_age1 <- subset(data2, data2$age<31)
data3_age1 <- subset(data3, data3$age<31)

#models by treatment type
lm1a <- lm_robust(comparison_1 ~ as.factor(trt) + education + income + newjob, data=data1_age1, se_type="HC3")
lm2a <- lm_robust(comparison_1 ~ as.factor(trt) + lr_1  + natid_covar1_1 + income + newjob, data=data2_age1, se_type="HC3")
lm3a <- lm_robust(comparison_1 ~ as.factor(trt) + lr_1 + unemployed, data=data3_age1, se_type="HC3")

idt <- c("Identity Threat\nvs. Control", lm1a$coefficients[2]) 
idt[2] <- lapply(idt[2], as.numeric); idt <- as.data.frame(idt)
colnames(idt) <- c("Treatment_Type", "Effect")
idt$se <- (lm1a$std.error)[2]
idt <- idt %>% mutate(cih90 = Effect + 1.645*se)
idt <- idt %>% mutate(cil90 = Effect - 1.645*se)
idt <- idt %>% mutate(cih95 = Effect + 1.96*se)
idt <- idt %>% mutate(cil95 = Effect - 1.96*se)

cit <- c("Whataboutism\nvs. Identity Threat", lm2a$coefficients[2]) 
cit[2] <- lapply(cit[2], as.numeric); cit <- as.data.frame(cit)
colnames(cit) <- c("Treatment_Type", "Effect")
cit$se <- (lm2a$std.error)[2]
cit <- cit %>% mutate(cih90 = Effect + 1.645*se)
cit <- cit %>% mutate(cil90 = Effect - 1.645*se)
cit <- cit %>% mutate(cih95 = Effect + 1.96*se)
cit <- cit %>% mutate(cil95 = Effect - 1.96*se)

com <- c("Whataboutism\nvs. Control", lm3a$coefficients[2]) 
com[2] <- lapply(com[2], as.numeric); com <- as.data.frame(com)
colnames(com) <- c("Treatment_Type", "Effect")
com$se <- (lm3a$std.error)[2]
com <- com %>% mutate(cih90 = Effect + 1.645*se)
com <- com %>% mutate(cil90 = Effect - 1.645*se)
com <- com %>% mutate(cih95 = Effect + 1.96*se)
com <- com %>% mutate(cil95 = Effect - 1.96*se)

#bind results by treatment type
results_age1 <- rbind(unlist(idt), unlist(cit), unlist(com))
colnames(results_age1) <- c("Treatment_Type", "Effect", "se", "cih90", "cil90", "cih95", "cil95")
results_age1 <- data.frame(results_age1)
results_age1$Effect <- as.numeric(results_age1$Effect)
results_age1$cil90 <- as.numeric(results_age1$cil90)
results_age1$cih90 <- as.numeric(results_age1$cih90)
results_age1$cil95 <- as.numeric(results_age1$cil95)
results_age1$cih95 <- as.numeric(results_age1$cih95)

results_age1$Age<- "18-30"

#age group 2

data1_age2 <- subset(data1, data1$age>30 & data1$age<41)
data2_age2 <- subset(data2, data2$age>30 & data2$age<41)
data3_age2 <- subset(data3, data3$age>30 & data3$age<41)

#models by treatment type
lm1a <- lm_robust(comparison_1 ~ as.factor(trt) + education + income + newjob, data=data1_age2, se_type="HC3")
lm2a <- lm_robust(comparison_1 ~ as.factor(trt) + lr_1  + natid_covar1_1 + income + newjob, data=data2_age2, se_type="HC3")
lm3a <- lm_robust(comparison_1 ~ as.factor(trt) + lr_1 + unemployed, data=data3_age2, se_type="HC3")

idt <- c("Identity Threat\nvs. Control", lm1a$coefficients[2]) 
idt[2] <- lapply(idt[2], as.numeric); idt <- as.data.frame(idt)
colnames(idt) <- c("Treatment_Type", "Effect")
idt$se <- (lm1a$std.error)[2]
idt <- idt %>% mutate(cih90 = Effect + 1.645*se)
idt <- idt %>% mutate(cil90 = Effect - 1.645*se)
idt <- idt %>% mutate(cih95 = Effect + 1.96*se)
idt <- idt %>% mutate(cil95 = Effect - 1.96*se)

cit <- c("Whataboutism\nvs. Identity Threat", lm2a$coefficients[2]) 
cit[2] <- lapply(cit[2], as.numeric); cit <- as.data.frame(cit)
colnames(cit) <- c("Treatment_Type", "Effect")
cit$se <- (lm2a$std.error)[2]
cit <- cit %>% mutate(cih90 = Effect + 1.645*se)
cit <- cit %>% mutate(cil90 = Effect - 1.645*se)
cit <- cit %>% mutate(cih95 = Effect + 1.96*se)
cit <- cit %>% mutate(cil95 = Effect - 1.96*se)

com <- c("Whataboutism\nvs. Control", lm3a$coefficients[2]) 
com[2] <- lapply(com[2], as.numeric); com <- as.data.frame(com)
colnames(com) <- c("Treatment_Type", "Effect")
com$se <- (lm3a$std.error)[2]
com <- com %>% mutate(cih90 = Effect + 1.645*se)
com <- com %>% mutate(cil90 = Effect - 1.645*se)
com <- com %>% mutate(cih95 = Effect + 1.96*se)
com <- com %>% mutate(cil95 = Effect - 1.96*se)

#bind results by treatment type
results_age2 <- rbind(unlist(idt), unlist(cit), unlist(com))
colnames(results_age2) <- c("Treatment_Type", "Effect", "se", "cih90", "cil90", "cih95", "cil95")
results_age2 <- data.frame(results_age2)
results_age2$Effect <- as.numeric(results_age2$Effect)
results_age2$cil90 <- as.numeric(results_age2$cil90)
results_age2$cih90 <- as.numeric(results_age2$cih90)
results_age2$cil95 <- as.numeric(results_age2$cil95)
results_age2$cih95 <- as.numeric(results_age2$cih95)
results_age2$Age<- "31-40"

#age group 3

data1_age3 <- subset(data1, data1$age>40 & data1$age<51)
data2_age3 <- subset(data2, data2$age>40 & data2$age<51)
data3_age3 <- subset(data3, data3$age>40 & data3$age<51)

#models by treatment type
lm1a <- lm_robust(comparison_1 ~ as.factor(trt) + education + income + newjob, data=data1_age3, se_type="HC3")
lm2a <- lm_robust(comparison_1 ~ as.factor(trt) + lr_1  + natid_covar1_1 + income + newjob, data=data2_age3, se_type="HC3")
lm3a <- lm_robust(comparison_1 ~ as.factor(trt) + lr_1 + unemployed, data=data3_age3, se_type="HC3")

idt <- c("Identity Threat\nvs. Control", lm1a$coefficients[2]) 
idt[2] <- lapply(idt[2], as.numeric); idt <- as.data.frame(idt)
colnames(idt) <- c("Treatment_Type", "Effect")
idt$se <- (lm1a$std.error)[2]
idt <- idt %>% mutate(cih90 = Effect + 1.645*se)
idt <- idt %>% mutate(cil90 = Effect - 1.645*se)
idt <- idt %>% mutate(cih95 = Effect + 1.96*se)
idt <- idt %>% mutate(cil95 = Effect - 1.96*se)

cit <- c("Whataboutism\nvs. Identity Threat", lm2a$coefficients[2]) 
cit[2] <- lapply(cit[2], as.numeric); cit <- as.data.frame(cit)
colnames(cit) <- c("Treatment_Type", "Effect")
cit$se <- (lm2a$std.error)[2]
cit <- cit %>% mutate(cih90 = Effect + 1.645*se)
cit <- cit %>% mutate(cil90 = Effect - 1.645*se)
cit <- cit %>% mutate(cih95 = Effect + 1.96*se)
cit <- cit %>% mutate(cil95 = Effect - 1.96*se)

com <- c("Whataboutism\nvs. Control", lm3a$coefficients[2]) 
com[2] <- lapply(com[2], as.numeric); com <- as.data.frame(com)
colnames(com) <- c("Treatment_Type", "Effect")
com$se <- (lm3a$std.error)[2]
com <- com %>% mutate(cih90 = Effect + 1.645*se)
com <- com %>% mutate(cil90 = Effect - 1.645*se)
com <- com %>% mutate(cih95 = Effect + 1.96*se)
com <- com %>% mutate(cil95 = Effect - 1.96*se)

#bind results by treatment type
results_age3 <- rbind(unlist(idt), unlist(cit), unlist(com))
colnames(results_age3) <- c("Treatment_Type", "Effect", "se", "cih90", "cil90", "cih95", "cil95")
results_age3 <- data.frame(results_age3)
results_age3$Effect <- as.numeric(results_age3$Effect)
results_age3$cil90 <- as.numeric(results_age3$cil90)
results_age3$cih90 <- as.numeric(results_age3$cih90)
results_age3$cil95 <- as.numeric(results_age3$cil95)
results_age3$cih95 <- as.numeric(results_age3$cih95)
results_age3$Age<- "41-50"

#age group 4

data1_age4 <- subset(data1, data1$age>50 & data1$age<61)
data2_age4 <- subset(data2, data2$age>50 & data2$age<61)
data3_age4 <- subset(data3, data3$age>50 & data3$age<61)

#models by treatment type
lm1a <- lm_robust(comparison_1 ~ as.factor(trt) + education + income + newjob, data=data1_age4, se_type="HC3")
lm2a <- lm_robust(comparison_1 ~ as.factor(trt) + lr_1  + natid_covar1_1 + income + newjob, data=data2_age4, se_type="HC3")
lm3a <- lm_robust(comparison_1 ~ as.factor(trt) + lr_1 + unemployed, data=data3_age4, se_type="HC3")

idt <- c("Identity Threat\nvs. Control", lm1a$coefficients[2]) 
idt[2] <- lapply(idt[2], as.numeric); idt <- as.data.frame(idt)
colnames(idt) <- c("Treatment_Type", "Effect")
idt$se <- (lm1a$std.error)[2]
idt <- idt %>% mutate(cih90 = Effect + 1.645*se)
idt <- idt %>% mutate(cil90 = Effect - 1.645*se)
idt <- idt %>% mutate(cih95 = Effect + 1.96*se)
idt <- idt %>% mutate(cil95 = Effect - 1.96*se)

cit <- c("Whataboutism\nvs. Identity Threat", lm2a$coefficients[2]) 
cit[2] <- lapply(cit[2], as.numeric); cit <- as.data.frame(cit)
colnames(cit) <- c("Treatment_Type", "Effect")
cit$se <- (lm2a$std.error)[2]
cit <- cit %>% mutate(cih90 = Effect + 1.645*se)
cit <- cit %>% mutate(cil90 = Effect - 1.645*se)
cit <- cit %>% mutate(cih95 = Effect + 1.96*se)
cit <- cit %>% mutate(cil95 = Effect - 1.96*se)

com <- c("Whataboutism\nvs. Control", lm3a$coefficients[2]) 
com[2] <- lapply(com[2], as.numeric); com <- as.data.frame(com)
colnames(com) <- c("Treatment_Type", "Effect")
com$se <- (lm3a$std.error)[2]
com <- com %>% mutate(cih90 = Effect + 1.645*se)
com <- com %>% mutate(cil90 = Effect - 1.645*se)
com <- com %>% mutate(cih95 = Effect + 1.96*se)
com <- com %>% mutate(cil95 = Effect - 1.96*se)

#bind results by treatment type
results_age4 <- rbind(unlist(idt), unlist(cit), unlist(com))
colnames(results_age4) <- c("Treatment_Type", "Effect", "se", "cih90", "cil90", "cih95", "cil95")
results_age4 <- data.frame(results_age4)
results_age4$Effect <- as.numeric(results_age4$Effect)
results_age4$cil90 <- as.numeric(results_age4$cil90)
results_age4$cih90 <- as.numeric(results_age4$cih90)
results_age4$cil95 <- as.numeric(results_age4$cil95)
results_age4$cih95 <- as.numeric(results_age4$cih95)
results_age4$Age<- "51-60"

#age group 5

data1_age5 <- subset(data1, data1$age>60)
data2_age5 <- subset(data2, data2$age>60)
data3_age5 <- subset(data3, data3$age>60)

#models by treatment type
lm1a <- lm_robust(comparison_1 ~ as.factor(trt) + education + income + newjob, data=data1_age5, se_type="HC3")
lm2a <- lm_robust(comparison_1 ~ as.factor(trt) + lr_1  + natid_covar1_1 + income + newjob, data=data2_age5, se_type="HC3")
lm3a <- lm_robust(comparison_1 ~ as.factor(trt) + lr_1 + unemployed, data=data3_age5, se_type="HC3")

idt <- c("Identity Threat\nvs. Control", lm1a$coefficients[2]) 
idt[2] <- lapply(idt[2], as.numeric); idt <- as.data.frame(idt)
colnames(idt) <- c("Treatment_Type", "Effect")
idt$se <- (lm1a$std.error)[2]
idt <- idt %>% mutate(cih90 = Effect + 1.645*se)
idt <- idt %>% mutate(cil90 = Effect - 1.645*se)
idt <- idt %>% mutate(cih95 = Effect + 1.96*se)
idt <- idt %>% mutate(cil95 = Effect - 1.96*se)

cit <- c("Whataboutism\nvs. Identity Threat", lm2a$coefficients[2]) 
cit[2] <- lapply(cit[2], as.numeric); cit <- as.data.frame(cit)
colnames(cit) <- c("Treatment_Type", "Effect")
cit$se <- (lm2a$std.error)[2]
cit <- cit %>% mutate(cih90 = Effect + 1.645*se)
cit <- cit %>% mutate(cil90 = Effect - 1.645*se)
cit <- cit %>% mutate(cih95 = Effect + 1.96*se)
cit <- cit %>% mutate(cil95 = Effect - 1.96*se)

com <- c("Whataboutism\nvs. Control", lm3a$coefficients[2]) 
com[2] <- lapply(com[2], as.numeric); com <- as.data.frame(com)
colnames(com) <- c("Treatment_Type", "Effect")
com$se <- (lm3a$std.error)[2]
com <- com %>% mutate(cih90 = Effect + 1.645*se)
com <- com %>% mutate(cil90 = Effect - 1.645*se)
com <- com %>% mutate(cih95 = Effect + 1.96*se)
com <- com %>% mutate(cil95 = Effect - 1.96*se)

#bind results by treatment type
results_age5 <- rbind(unlist(idt), unlist(cit), unlist(com))
colnames(results_age5) <- c("Treatment_Type", "Effect", "se", "cih90", "cil90", "cih95", "cil95")
results_age5 <- data.frame(results_age5)
results_age5$Effect <- as.numeric(results_age5$Effect)
results_age5$cil90 <- as.numeric(results_age5$cil90)
results_age5$cih90 <- as.numeric(results_age5$cih90)
results_age5$cil95 <- as.numeric(results_age5$cil95)
results_age5$cih95 <- as.numeric(results_age5$cih95)
results_age5$Age<- "61+"

#combine age groups and plot

age_hte <- rbind(results_age1, results_age2, results_age3, results_age4, results_age5)
age_hte$Age <- factor(age_hte$Age, levels = c("18-30", "31-40", "41-50", "51-60", "61+"))

wb2 <- ggplot(age_hte, aes(x = Age, y = Effect, colour=factor(Treatment_Type, level=factor_levels))) +
  geom_point(stat="identity", position=position_dodge(width=0.5), size=3) + 
  labs(tag = "B") +
  scale_y_continuous(limits = c(-1.55, 1.55)) +
  geom_linerange(aes(ymin = cil95, ymax = cih95), 
                 position = position_dodge(width=0.5,),
                 size = 1) +
  geom_linerange(aes(ymin = cil90, ymax = cih90), 
                 position = position_dodge(width=0.5),
                 size = 2) +
  ylab("Treatment Effect") +
  xlab("Age") +
  geom_hline(yintercept=0, linetype="dashed", color="blue") +
  labs(color='Update Type') +
  theme_bw()+
  theme(
    axis.text=element_text(size=16),
    plot.title = element_text(size=25),
    axis.title=element_text(size=18),
    legend.text=element_text(size=12),
    legend.title=element_text(size=16),
    legend.key.height = unit(1.5, 'cm')
  )

####TABLE D6(C) - Whataboutism by Gender####

# women
data1_f <- subset(data1, data1$female == 1)
data2_f <- subset(data2, data2$female == 1)
data3_f <- subset(data3, data3$female == 1)

#models by treatment type
lm1a <- lm_robust(comparison_1 ~ as.factor(trt) + education + income + newjob, data=data1_f, se_type="HC3")
lm2a <- lm_robust(comparison_1 ~ as.factor(trt) + lr_1 + natid_covar1_1 + income + newjob, data = data2_f, se_type="HC3")
lm3a <- lm_robust(comparison_1 ~ as.factor(trt) + lr_1 + unemployed, data = data3_f, se_type="HC3")

idt <- data.frame(
  Treatment_Type = "Identity Threat\nvs. Control",
  Effect = lm1a$coefficients[2],
  se = lm1a$std.error[2]
)
idt <- idt %>% 
  mutate(
    ci90h = Effect + 1.645 * se,
    ci90l = Effect - 1.645 * se,
    ci95h = Effect + 1.96 * se,
    ci95l = Effect - 1.96 * se
  )

cit <- data.frame(
  Treatment_Type = "Whataboutism\nvs. Identity Threat",
  Effect = lm2a$coefficients[2],
  se = lm2a$std.error[2]
)
cit <- cit %>% 
  mutate(
    ci90h = Effect + 1.645 * se,
    ci90l = Effect - 1.645 * se,
    ci95h = Effect + 1.96 * se,
    ci95l = Effect - 1.96 * se
  )

com <- data.frame(
  Treatment_Type = "Whataboutism\nvs. Control",
  Effect = lm3a$coefficients[2],
  se = lm3a$std.error[2]
)
com <- com %>% 
  mutate(
    ci90h = Effect + 1.645 * se,
    ci90l = Effect - 1.645 * se,
    ci95h = Effect + 1.96 * se,
    ci95l = Effect - 1.96 * se
  )

#bind results by treatment type
results_f <- rbind(idt, cit, com)
results_f$Gender <- "Female"

# men
data1_m <- subset(data1, data1$female == 0)
data2_m <- subset(data2, data2$female == 0)
data3_m <- subset(data3, data3$female == 0)

#models by treatment type
lm1a <- lm_robust(comparison_1 ~ as.factor(trt) + education + income + newjob, data=data1_m, se_type="HC3")
lm2a <- lm_robust(comparison_1 ~ as.factor(trt) + lr_1 + natid_covar1_1 + income + newjob, data = data2_m, se_type="HC3")
lm3a <- lm_robust(comparison_1 ~ as.factor(trt) + lr_1 + unemployed, data = data3_m, se_type="HC3")

idt <- data.frame(
  Treatment_Type = "Identity Threat\nvs. Control",
  Effect = lm1a$coefficients[2],
  se = lm1a$std.error[2]
)
idt <- idt %>% 
  mutate(
    ci90h = Effect + 1.645 * se,
    ci90l = Effect - 1.645 * se,
    ci95h = Effect + 1.96 * se,
    ci95l = Effect - 1.96 * se
  )

cit <- data.frame(
  Treatment_Type = "Whataboutism\nvs. Identity Threat",
  Effect = lm2a$coefficients[2],
  se = lm2a$std.error[2]
)
cit <- cit %>% 
  mutate(
    ci90h = Effect + 1.645 * se,
    ci90l = Effect - 1.645 * se,
    ci95h = Effect + 1.96 * se,
    ci95l = Effect - 1.96 * se
  )

com <- data.frame(
  Treatment_Type = "Whataboutism\nvs. Control",
  Effect = lm3a$coefficients[2],
  se = lm3a$std.error[2]
)
com <- com %>% 
  mutate(
    ci90h = Effect + 1.645 * se,
    ci90l = Effect - 1.645 * se,
    ci95h = Effect + 1.96 * se,
    ci95l = Effect - 1.96 * se
  )

#bind results by treatment type
results_m <- rbind(idt, cit, com)
results_m$Gender <- "Male"

# Combine women/men and plot
gender_hte <- rbind(results_f, results_m)
gender_hte$Gender <- factor(gender_hte$Gender, levels = c("Female", "Male"))

wb3 <- ggplot(gender_hte, aes(x = Gender, y = Effect, colour=factor(Treatment_Type, level=factor_levels))) +
  geom_point(stat="identity", position=position_dodge(width=0.3), size=3) + 
  labs(tag = "C") +
  scale_y_continuous(limits = c(-1.55, 1.55)) +
  geom_linerange(aes(ymin = ci95l, ymax = ci95h), 
                 position = position_dodge(width=0.3,),
                 size = 1) +
  geom_linerange(aes(ymin = ci90l, ymax = ci90h), 
                 position = position_dodge(width=0.3),
                 size = 2)+
  ylab("Treatment Effect") +
  xlab("Gender") +
  geom_hline(yintercept=0, linetype="dashed", color="blue") +
  labs(color='Update Type') +
  theme_bw()+
  theme(
    axis.text=element_text(size=16),
    plot.title = element_text(size=25),
    axis.title=element_text(size=18),
    legend.text=element_text(size=12),
    legend.title=element_text(size=16),
    legend.key.height = unit(1.5, 'cm')
  )


####TABLE D6(D) - Whataboutism by Pre-Treatment National Identification####

#quintile 1

data1_nid1 <- subset(data1, data1$nat_quint==1)
data2_nid1 <- subset(data2, data2$nat_quint==1)
data3_nid1 <- subset(data3, data3$nat_quint==1)

#models by treatment type
lm1a <- lm_robust(comparison_1 ~ as.factor(trt) + education + income + newjob, data=data1_nid1, se_type="HC3")
lm2a <- lm_robust(comparison_1 ~ as.factor(trt) + lr_1 + natid_covar1_1 + income + newjob, data = data2_nid1, se_type="HC3")
lm3a <- lm_robust(comparison_1 ~ as.factor(trt) + lr_1 + unemployed, data = data3_nid1, se_type="HC3")

idt <- c("Identity Threat\nvs. Control", lm1a$coefficients[2]) 
idt[2] <- lapply(idt[2], as.numeric); idt <- as.data.frame(idt)
colnames(idt) <- c("Treatment_Type", "Effect")
idt$se <- (lm1a$std.error)[2]
idt <- idt %>% mutate(cih95 = Effect + 1.96*se)
idt <- idt %>% mutate(cil95 = Effect - 1.96*se)
idt <- idt %>% mutate(cih90 = Effect + 1.645*se)
idt <- idt %>% mutate(cil90 = Effect - 1.645*se)

cit <- c("Whataboutism\nvs. Identity Threat", lm2a$coefficients[2]) 
cit[2] <- lapply(cit[2], as.numeric); cit <- as.data.frame(cit)
colnames(cit) <- c("Treatment_Type", "Effect")
cit$se <- (lm2a$std.error)[2]
cit <- cit %>% mutate(cih95 = Effect + 1.96*se)
cit <- cit %>% mutate(cil95 = Effect - 1.96*se)
cit <- cit %>% mutate(cih90 = Effect + 1.645*se)
cit <- cit %>% mutate(cil90 = Effect - 1.645*se)

com <- c("Whataboutism\nvs. Control", lm3a$coefficients[2]) 
com[2] <- lapply(com[2], as.numeric); com <- as.data.frame(com)
colnames(com) <- c("Treatment_Type", "Effect")
com$se <- (lm3a$std.error)[2]
com <- com %>% mutate(cih95 = Effect + 1.96*se)
com <- com %>% mutate(cil95 = Effect - 1.96*se)
com <- com %>% mutate(cih90 = Effect + 1.645*se)
com <- com %>% mutate(cil90 = Effect - 1.645*se)

#bind results by treatment type
results_nid1 <- rbind(unlist(idt), unlist(cit), unlist(com))
colnames(results_nid1) <- c("Treatment_Type", "Effect", "se", "cih95", "cil95", "cih90", "cil90")
results_nid1 <- data.frame(results_nid1)
results_nid1$Effect <- as.numeric(results_nid1$Effect)
results_nid1$cil95 <- as.numeric(results_nid1$cil95)
results_nid1$cil90 <- as.numeric(results_nid1$cil90)
results_nid1$cih95 <- as.numeric(results_nid1$cih95)
results_nid1$cih90 <- as.numeric(results_nid1$cih90)
results_nid1$Quintile<- 1

#quintile 2

data1_nid2 <- subset(data1, data1$nat_quint==2)
data2_nid2 <- subset(data2, data2$nat_quint==2)
data3_nid2 <- subset(data3, data3$nat_quint==2)

#models by treatment type
lm1a <- lm_robust(comparison_1 ~ as.factor(trt) + education + income + newjob, data=data1_nid2, se_type="HC3")
lm2a <- lm_robust(comparison_1 ~ as.factor(trt) + lr_1 + natid_covar1_1 + income + newjob, data = data2_nid2, se_type="HC3")
lm3a <- lm_robust(comparison_1 ~ as.factor(trt) + lr_1 + unemployed, data = data3_nid2, se_type="HC3")

idt <- c("Identity Threat\nvs. Control", lm1a$coefficients[2]) 
idt[2] <- lapply(idt[2], as.numeric); idt <- as.data.frame(idt)
colnames(idt) <- c("Treatment_Type", "Effect")
idt$se <- (lm1a$std.error)[2]
idt <- idt %>% mutate(cih95 = Effect + 1.96*se)
idt <- idt %>% mutate(cil95 = Effect - 1.96*se)
idt <- idt %>% mutate(cih90 = Effect + 1.645*se)
idt <- idt %>% mutate(cil90 = Effect - 1.645*se)

cit <- c("Whataboutism\nvs. Identity Threat", lm2a$coefficients[2]) 
cit[2] <- lapply(cit[2], as.numeric); cit <- as.data.frame(cit)
colnames(cit) <- c("Treatment_Type", "Effect")
cit$se <- (lm2a$std.error)[2]
cit <- cit %>% mutate(cih95 = Effect + 1.96*se)
cit <- cit %>% mutate(cil95 = Effect - 1.96*se)
cit <- cit %>% mutate(cih90 = Effect + 1.645*se)
cit <- cit %>% mutate(cil90 = Effect - 1.645*se)

com <- c("Whataboutism\nvs. Control", lm3a$coefficients[2]) 
com[2] <- lapply(com[2], as.numeric); com <- as.data.frame(com)
colnames(com) <- c("Treatment_Type", "Effect")
com$se <- (lm3a$std.error)[2]
com <- com %>% mutate(cih95 = Effect + 1.96*se)
com <- com %>% mutate(cil95 = Effect - 1.96*se)
com <- com %>% mutate(cih90 = Effect + 1.645*se)
com <- com %>% mutate(cil90 = Effect - 1.645*se)

#bind results by treatment type
results_nid2 <- rbind(unlist(idt), unlist(cit), unlist(com))
colnames(results_nid2) <- c("Treatment_Type", "Effect", "se", "cih95", "cil95", "cih90", "cil90")
results_nid2 <- data.frame(results_nid2)
results_nid2$Effect <- as.numeric(results_nid2$Effect)
results_nid2$cil95 <- as.numeric(results_nid2$cil95)
results_nid2$cil90 <- as.numeric(results_nid2$cil90)
results_nid2$cih95 <- as.numeric(results_nid2$cih95)
results_nid2$cih90 <- as.numeric(results_nid2$cih90)
results_nid2$Quintile<- 2

#quintile 3

data1_nid3 <- subset(data1, data1$nat_quint==3)
data2_nid3 <- subset(data2, data2$nat_quint==3)
data3_nid3 <- subset(data3, data3$nat_quint==3)

#models by treatment type
lm1a <- lm_robust(comparison_1 ~ as.factor(trt) + education + income + newjob, data=data1_nid3, se_type="HC3")
lm2a <- lm_robust(comparison_1 ~ as.factor(trt) + lr_1 + natid_covar1_1 + income + newjob, data = data2_nid3, se_type="HC3")
lm3a <- lm_robust(comparison_1 ~ as.factor(trt) + lr_1 + unemployed, data = data3_nid3, se_type="HC3")

idt <- c("Identity Threat\nvs. Control", lm1a$coefficients[2]) 
idt[2] <- lapply(idt[2], as.numeric); idt <- as.data.frame(idt)
colnames(idt) <- c("Treatment_Type", "Effect")
idt$se <- (lm1a$std.error)[2]
idt <- idt %>% mutate(cih95 = Effect + 1.96*se)
idt <- idt %>% mutate(cil95 = Effect - 1.96*se)
idt <- idt %>% mutate(cih90 = Effect + 1.645*se)
idt <- idt %>% mutate(cil90 = Effect - 1.645*se)

cit <- c("Whataboutism\nvs. Identity Threat", lm2a$coefficients[2]) 
cit[2] <- lapply(cit[2], as.numeric); cit <- as.data.frame(cit)
colnames(cit) <- c("Treatment_Type", "Effect")
cit$se <- (lm2a$std.error)[2]
cit <- cit %>% mutate(cih95 = Effect + 1.96*se)
cit <- cit %>% mutate(cil95 = Effect - 1.96*se)
cit <- cit %>% mutate(cih90 = Effect + 1.645*se)
cit <- cit %>% mutate(cil90 = Effect - 1.645*se)

com <- c("Whataboutism\nvs. Control", lm3a$coefficients[2]) 
com[2] <- lapply(com[2], as.numeric); com <- as.data.frame(com)
colnames(com) <- c("Treatment_Type", "Effect")
com$se <- (lm3a$std.error)[2]
com <- com %>% mutate(cih95 = Effect + 1.96*se)
com <- com %>% mutate(cil95 = Effect - 1.96*se)
com <- com %>% mutate(cih90 = Effect + 1.645*se)
com <- com %>% mutate(cil90 = Effect - 1.645*se)

#bind results by treatment type
results_nid3 <- rbind(unlist(idt), unlist(cit), unlist(com))
colnames(results_nid3) <- c("Treatment_Type", "Effect", "se", "cih95", "cil95", "cih90", "cil90")
results_nid3 <- data.frame(results_nid3)
results_nid3$Effect <- as.numeric(results_nid3$Effect)
results_nid3$cil95 <- as.numeric(results_nid3$cil95)
results_nid3$cil90 <- as.numeric(results_nid3$cil90)
results_nid3$cih95 <- as.numeric(results_nid3$cih95)
results_nid3$cih90 <- as.numeric(results_nid3$cih90)
results_nid3$Quintile<- 3

#quintile 4

data1_nid4 <- subset(data1, data1$nat_quint==4)
data2_nid4 <- subset(data2, data2$nat_quint==4)
data3_nid4 <- subset(data3, data3$nat_quint==4)

#models by treatment type
lm1a <- lm_robust(comparison_1 ~ as.factor(trt) + education + income + newjob, data=data1_nid4, se_type="HC3")
lm2a <- lm_robust(comparison_1 ~ as.factor(trt) + lr_1 + natid_covar1_1 + income + newjob, data = data2_nid4, se_type="HC3")
lm3a <- lm_robust(comparison_1 ~ as.factor(trt) + lr_1 + unemployed, data = data3_nid4, se_type="HC3")

idt <- c("Identity Threat\nvs. Control", lm1a$coefficients[2]) 
idt[2] <- lapply(idt[2], as.numeric); idt <- as.data.frame(idt)
colnames(idt) <- c("Treatment_Type", "Effect")
idt$se <- (lm1a$std.error)[2]
idt <- idt %>% mutate(cih95 = Effect + 1.96*se)
idt <- idt %>% mutate(cil95 = Effect - 1.96*se)
idt <- idt %>% mutate(cih90 = Effect + 1.645*se)
idt <- idt %>% mutate(cil90 = Effect - 1.645*se)

cit <- c("Whataboutism\nvs. Identity Threat", lm2a$coefficients[2]) 
cit[2] <- lapply(cit[2], as.numeric); cit <- as.data.frame(cit)
colnames(cit) <- c("Treatment_Type", "Effect")
cit$se <- (lm2a$std.error)[2]
cit <- cit %>% mutate(cih95 = Effect + 1.96*se)
cit <- cit %>% mutate(cil95 = Effect - 1.96*se)
cit <- cit %>% mutate(cih90 = Effect + 1.645*se)
cit <- cit %>% mutate(cil90 = Effect - 1.645*se)

com <- c("Whataboutism\nvs. Control", lm3a$coefficients[2]) 
com[2] <- lapply(com[2], as.numeric); com <- as.data.frame(com)
colnames(com) <- c("Treatment_Type", "Effect")
com$se <- (lm3a$std.error)[2]
com <- com %>% mutate(cih95 = Effect + 1.96*se)
com <- com %>% mutate(cil95 = Effect - 1.96*se)
com <- com %>% mutate(cih90 = Effect + 1.645*se)
com <- com %>% mutate(cil90 = Effect - 1.645*se)

#bind results by treatment type
results_nid4 <- rbind(unlist(idt), unlist(cit), unlist(com))
colnames(results_nid4) <- c("Treatment_Type", "Effect", "se", "cih95", "cil95", "cih90", "cil90")
results_nid4 <- data.frame(results_nid4)
results_nid4$Effect <- as.numeric(results_nid4$Effect)
results_nid4$cil95 <- as.numeric(results_nid4$cil95)
results_nid4$cil90 <- as.numeric(results_nid4$cil90)
results_nid4$cih95 <- as.numeric(results_nid4$cih95)
results_nid4$cih90 <- as.numeric(results_nid4$cih90)
results_nid4$Quintile<- 4

#quintile 5

data1_nid5 <- subset(data1, data1$nat_quint==5)
data2_nid5 <- subset(data2, data2$nat_quint==5)
data3_nid5 <- subset(data3, data3$nat_quint==5)

#models by treatment type
lm1a <- lm_robust(comparison_1 ~ as.factor(trt) + education + income + newjob, data=data1_nid5, se_type="HC3")
lm2a <- lm_robust(comparison_1 ~ as.factor(trt) + lr_1 + natid_covar1_1 + income + newjob, data = data2_nid5, se_type="HC3")
lm3a <- lm_robust(comparison_1 ~ as.factor(trt) + lr_1 + unemployed, data = data3_nid5, se_type="HC3")

idt <- c("Identity Threat\nvs. Control", lm1a$coefficients[2]) 
idt[2] <- lapply(idt[2], as.numeric); idt <- as.data.frame(idt)
colnames(idt) <- c("Treatment_Type", "Effect")
idt$se <- (lm1a$std.error)[2]
idt <- idt %>% mutate(cih95 = Effect + 1.96*se)
idt <- idt %>% mutate(cil95 = Effect - 1.96*se)
idt <- idt %>% mutate(cih90 = Effect + 1.645*se)
idt <- idt %>% mutate(cil90 = Effect - 1.645*se)

cit <- c("Whataboutism\nvs. Identity Threat", lm2a$coefficients[2]) 
cit[2] <- lapply(cit[2], as.numeric); cit <- as.data.frame(cit)
colnames(cit) <- c("Treatment_Type", "Effect")
cit$se <- (lm2a$std.error)[2]
cit <- cit %>% mutate(cih95 = Effect + 1.96*se)
cit <- cit %>% mutate(cil95 = Effect - 1.96*se)
cit <- cit %>% mutate(cih90 = Effect + 1.645*se)
cit <- cit %>% mutate(cil90 = Effect - 1.645*se)

com <- c("Whataboutism\nvs. Control", lm3a$coefficients[2]) 
com[2] <- lapply(com[2], as.numeric); com <- as.data.frame(com)
colnames(com) <- c("Treatment_Type", "Effect")
com$se <- (lm3a$std.error)[2]
com <- com %>% mutate(cih95 = Effect + 1.96*se)
com <- com %>% mutate(cil95 = Effect - 1.96*se)
com <- com %>% mutate(cih90 = Effect + 1.645*se)
com <- com %>% mutate(cil90 = Effect - 1.645*se)

#bind results by treatment type
results_nid5 <- rbind(unlist(idt), unlist(cit), unlist(com))
colnames(results_nid5) <- c("Treatment_Type", "Effect", "se", "cih95", "cil95", "cih90", "cil90")
results_nid5 <- data.frame(results_nid5)
results_nid5$Effect <- as.numeric(results_nid5$Effect)
results_nid5$cil95 <- as.numeric(results_nid5$cil95)
results_nid5$cil90 <- as.numeric(results_nid5$cil90)
results_nid5$cih95 <- as.numeric(results_nid5$cih95)
results_nid5$cih90 <- as.numeric(results_nid5$cih90)
results_nid5$Quintile<- 5

#combine quintiles and plot
nid_hte <- rbind(results_nid1, results_nid2, results_nid3, results_nid4, results_nid5)
nid_hte$Effect <- as.numeric(as.character(nid_hte$Effect))

wb4 <- ggplot(nid_hte, aes(x = Quintile, y = Effect, colour=factor(Treatment_Type, level=factor_levels))) +
  geom_point(stat="identity", position=position_dodge(width=0.5), size=3) + 
  labs(tag = "D") +
  scale_y_continuous(limits = c(-1.55, 1.55)) +
  geom_linerange(aes(ymin = cil95, ymax = cih95), 
                 position = position_dodge(width=0.5,),
                 size = 1) +
  geom_linerange(aes(ymin = cil90, ymax = cih90), 
                 position = position_dodge(width=0.5),
                 size = 2)+
  ylab("Treatment Effect") +
  xlab("National Identification") +
  geom_hline(yintercept=0, linetype="dashed", color="blue") +
  labs(color='Update Type') +
  theme_bw()+
  theme(
    axis.text=element_text(size=16),
    plot.title = element_text(size=25),
    axis.title=element_text(size=18),
    legend.text=element_text(size=12),
    legend.title=element_text(size=16),
    legend.key.height = unit(1.5, 'cm')
  )

####Combined Panel for Figure D6####

(wb1 | wb2) /
  (wb3 | wb4)+
  plot_layout(guides = "collect") &
  theme(legend.position = "right")

ggsave("figD6.png", width = 12, height = 10, dpi = 300)